/* News shock puzzle analysis */
/* By Toshihiro OKADA */


//** the time period of the simulated siries: ts=start te=end **
ts=1001; te=1190;
inds=5501; inde=5691;


//***** Choosing a set of parameters (posterior modes): Basic NK, Standard NK, and EPD-NK models *****//  
@#define data_techzero=1 
//* 0: Basic or Standard NK                      *//
//* 1: EPD-NK                                    *// 

@# define sum_shocks=1
//** NOt used **//

//********** Model choice ********* //
@#define tech=1
//* 1: EPD NK                *// 
//* 0: Basic or Standard NK  *//

//**** Shutting down endogenous technology for EPD-NK model ***// 
@#define ashock=1
//* 1: with endogenous tech.     *//
//* 0: without endogenous tech.  *//
//* 2: arbitrary value for "alfa" (R&D elasticity)  *//

//***** Price index *****//
@#define pindex=0
//* pindex==0: no price index  *//
//* pindex==1: with price index (use this for standard NK model with tech=0 and tech=0 *//
//*******************************//

//***** Time to Innovate (TTB) *****//
@#define TI=1
//* TI=1: 8-period Time to Innovate *//
//* TI=16: 16-period TTB  *//
//* TI=0: 1-period TTB, i.e., No Time to Innovate *//
//* Note when 16- period TTB is chosen, need to increase "phi" very slightly to satisfied the BK condion and chose phiaj_=1.035 in the some parameter adjustemnt section below **// 

//***** TTB period specification *****//
@#define ttb_spec=8
//**ttb_spec =8: 8-period time to innovate, ttb_spec=16: 16-period time to innovate **//
//* Note when 16- period TTB is chosen, need to increase "phi" very slightly to satisfied the BK condion and chose phiaj_=1.035 in the some parameter adjustemnt section below **// 

//********** Some parameter adjustment ********//
tfpwt_ =1.0;
alfaj_=1.0; 
phiaj_=1.0;  
psiaj_=1; 
rdscale_=1; 
rd_rt_=1; 

//***** Calvo Pricing Effect *****//
@#define calvo_pra=0
rhoaj_=0.5;  //* do not change this: setting rho (Calvo princing) to half of the baseline value *//
//* calvo_pra=0: baseline (posterior mode)  *//
//* calvo_pra=1: half of the baseline value *//

@#define eshock=16
//** Not used **//

//***** Shock Specification *****//
@#define ushock=1
@#define news=1
@#define dshock=0
//** dshock is not used **//
@#define cshock=1
@#define hshock=1
@#define ishock=1
@#define fshock=1
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//* ushock: ushock=0 shut down technology(u)shocks, 2: news: news=0, shut down  news shock but leave u shocks, cshock: consumption shock, hshock: hour shocks, ishock:interest rate shock *// 
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


//** croos correaltion of shocks **//
@#define crosscrr=1 


//**IRF figures spec: choose "0" for a usual use **//
@#define monte_bays_irf=0
@#define manu_irf=1

//**IRF CI bands for VAR: choose "1" for 90 % CI and "2" for 1 std bands **//
@#define irf_bands=2

@#define news_period=16
/* Not used */ 

@#define bayesian_irf=0
/* Not used  */

@#define bays_irf=0
/** This is for THEORETICAL multiple news shock IRFs **/


@#define bays_irf_non_news=0
/* Not used */


@#define check=0
/** Not used **/

@#define monte_irf=0
/** Not used **/

@#define sim_var=1
/** Not used: alsway choose 1 **/


//** Unrelaized shocks simulation **//
@#define unrealized=1
//** Not used **// 

@#define unrealized_news_sim=0 
//** Not used **// 

//*****************//
//****VAR METHOD**//
//****************//
@#define var_est_method=0
/** VAR on simultaed data: 2=bootstrap,  1=bayesian, 0=OLS or single bayesian **/
/** Not used: always choose 0 **/ 

@#define var_shock=0
 /** Not used **/

@#define median_mean=2
/** 1=median, 2=mean **/

@#define var_data=22
/** Not used **/

@#define conf_var=2
//** Not used **//


@#define sim_unrealized=0
//** Not used DO NOT TOUCH THIS(choose 0 or 2)**//

@#define method=1
//* NOT used *//

@#define mx_unrealize=1
/* mx_unrealized=1: some of expectations are unrealized (if wt=2 is chosen, 50 % of expectations are unrelalized, */
/* if wt=1 is chosen, 100 % are unrealized),  mx_unrealized=0: totally unrealized */
  
wt_=0.8;
//* Not used *//

//*******Not used option***//
@#define bps=0
/* Bandpass filtering: bps=1: band pass filtering is "on" */
@#define cmn_trd=0
/* common trend for Bpass */
@#define ctrd_spec=1
/* common trend specification: 1=no invest, 2=with invest */
@#define hp=0
@#define trend_sigma=0
//***********************************//

////** Line specification for figures **////////////
lsp = 1; /* setting a line specification for figures: 1='--', 2='-.', 3=':' and 4='-' */

/////** Figure Number **/////
@#define fignum =500
//* choose 400 / 500 *////



@#define tfp_errors=1
//** Not used**//

////*** Capital , habit & inv adj cost *****//  
@#define simple=0
//* =0: With capital, habit and inv adj cost *//
//* =1: No capital, no habit and no inv adj cost *//
//* =2: With capital & no habit and no inv adj cost *//
//* =3: With capital & inv adj cost, but No habit *//
//* =4: With capital & habit, NO inv adj cost *//
//***************************************//


@#define growth=0
//**Not used** //

@#define std_spec =0
//* std_spec=1: arbitrary std*//
//* std_spec=0: estimated std *//
//** Not used **//

//******** Conditional varinace; on-off ****//
@#define cndvar=0
/** cndvar=1: for (theoretical) conditional variance calculation (std=estimated values) **/  
/** cndvar=0: for otherwisee **/
//*******************************************//

@#define fig =0
/** fig=1: for thoeoretical simulation **/
/** fig=0: for monte carlo simulations **/  
//** Not used**/  

@#define fig_series=0
@#define full_fig=0




//**********************//
//** Loading data *********//
//***************************//

@#if data_techzero==1
  
 @#if eshock==16
 load('subdraws_phi60stv90mm5_mode.mat', 'aa0');
 @#endif
 
@#endif

@#if data_techzero==0
 
  @#if eshock==16

   //* Without Indexization *// 
   @#if pindex==0
   load('subdraws_BNKphi60stv90m3_mode.mat', 'aa0');
   @#endif

  //* With Indexization *// 
   @#if pindex==1
   load('subdraws_BNKphi60stv90index_mode.mat', 'aa0');
   @#endif
  @#endif

@#endif
//***********************************************//

//**** CHOICE OF SAMPLE PERIODS****//
@#define large_sample=0
sim_period=190;
//********************************//

me_std=0.00; 
mee_std=0.0; 
meee_std=0.0;

//**********************************//
//********SETTING SEEDS***********//
//**********************************//

//**** Monte carlo multiple news shocks IRF Simulation: random shock seeds *****//
//**** For multiple monte carlo news shock choose, e.g,dma=1 and dmb=500 *****//
//**** For multiple monte carlo news shock without bayesian coeff, choose ja=1001 and jb=1001, i.e., shut down bayesian random coeff. *****//
dma=1;
dmb=1;
//**********************************************************************//


//**** Random shock seeds for simulations of var, corr and etc: Use also for monte-carlo VAR simu (for sim_var=1)*****//
iia=1; /* start of random shock seed # */
iib=3000; /* end of random shock seed # */
//**************************************************************//

//**** Not used : Bayesian coefficient estimates random shock seeds ***********************//
ja=1;/* start of sub draws random shock seed #: FOR THE POSTERIOR Mode, CHOOSE 1001 */
jb=2;/* end of sub draws random shock seed # FOR THE POSTERIOR Mode, CHOOSE 1001 */
//**************************************************************//

//**** Not used: random shock seeds *****// 
jrdma=1; /* seed for news shock simulation random numbering*/
jrdmb=1; /* seed for news shock simulation random numbering*/
//**********************************************************************//

alfa_n = 150; /* to shut down technology effect */



@#define fullshock=1
//** Not used **//


ex_time_=4; 
//* Not used  *// 
@#define ex_time_periods=16
//* Not used  *//

@#if fullshock==1
 @#define shki =1
 @#define shku =1
 @#define shkuc =1
 @#define shkuh =1
 @#define shkue =1
 @#define shkf =1

 @#define shkue1 =1
 @#define shkue2 =1
 @#define shkue3 =1
 @#define shkue4 =1
 @#define shkue5 =1
 @#define shkue6 =1
 @#define shkue7 =1
 @#define shkue8 =1
 @#define shkue9 =1 
 @#if eshock==16
   @#define shkue10 =1 
   @#define shkue11 =1 
   @#define shkue12 =1 
   @#define shkue13 =1 
   @#define shkue14 =1 
   @#define shkue15 =1
  @#define shkue =1
 @#define shkuc =1
 @#define shkuh =1
@#define shkf =1
  @#endif  

 @#endif 


 @#define shock_nm=1
//** not used**//



for j=ja:jb

@#if tech==1
  @#if var_data==20
 var y_g, c_g, tiv_g, iv_g, h_g, w_g, tfp_g, gg, ggn, uf, pai,y,tfp,tfpe,tfpu,u,a, conf, conf12, conf11, conf10, conf9, conf8,c,k,iv,h,w,r,rr,p,ppai,m,pf,pf_exp,i,lamda,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,pfn,in,lamdan,yg,hg,fpai,mcgap,rd,rdn,an,nf,nfn,pinf,uc,um,uh,tiv, tivn ynobs,mcgape, mc, mcn  mce, fpaie, paie, paiee,ppaie, inftr,  tfpn, tq,tqn, shu, shi, shc, shh, shue1 shue3 shue4 shue8  aa aan sp stk  ;
 @#endif
 @#if var_data==21
 var y_g, c_g, tiv_g,iv_g,  h_g, w_g, tfp_g, gg, ggn, uf, pai,y,tfp,tfpu,u,a, conf, conf12, conf11, conf10, conf9, conf8,c,k,iv,h,w,r,rr,p,ppai,m,pf,pf_exp,i,lamda,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,pfn,in,lamdan,yg,hg,fpai,mcgap,rd,rdn,an,nf,nfn,pinf,uc,um,uh,tiv, tivn ynobs,mcgape, mc, mcn  mce, fpaie, paie, paiee,ppaie, inftr,  tfpn, tq,tqn, shu, shi, shc, shh, shue1 shue3 shue4 shue8  aa aan sp stk  ;
 @#endif
 @#if var_data==22
 var y_g, c_g, tiv_g,iv_g,  h_g, w_g, tfp_g, gg, ggn, uf, pai,y,tfp,tfpe,tfpu,u,a, conf, conf12, conf11, conf10, conf9, conf8,c,k,iv,h,w,r,rr,p,ppai,m,pf,pf_exp,i,lamda,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,pfn,in,lamdan,yg,hg,fpai,mcgap,rd,rdn,an,nf,nfn,pinf,uc,um,uh,tiv, tivn ynobs,mcgape, mc, mcn  mce, fpaie, paie, paiee,ppaie, inftr,  tfpn, tq,tqn, shu, shi, shc, shh, shue1 shue3 shue4 shue8  aa aan sp stk  ;
 @#endif
 @#if var_data==25
 var y_g, c_g, tiv_g, iv_g, h_g, w_g, tfp_g, gg, ggn, uf, pai,y,tfp,tfpu,u,a, conf, conf12, conf11, conf10, conf9, conf8,c,k,iv,h,w,r,rr,p,ppai,m,pf,pf_exp,i,lamda,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,pfn,in,lamdan,yg,hg,fpai,mcgap,rd,rdn,an,nf,nfn,pinf,uc,um,uh,tiv, tivn ynobs,mcgape, mc, mcn  mce, fpaie, paie, paiee,ppaie, inftr,  tfpn, tq,tqn, shu, shi, shc, shh, shue1 shue3 shue4 shue8  aa aan sp stk  ;
 @#endif
 

@#if eshock==16
  varexo ei eu euc euh ef eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11 eue12 eue13 eue14 eue15 eue etfp me mee meee;
@#endif 

@#if ttb_spec==8
parameters rdscale, tfpwt,tau,  hb,  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, psi, epsilon,   alfa, rhoum,rhouc, rhouh, rhou,rhoud,rhouf, sigc, sigh, sigm, ss, rho_eta, eta1, eta2, eta3, eta4,eta5,eta6,eta7,eta8,iss, rss, yss, kss, css,wss,mss, rdss, pfss,  hss, ivss, lamdass, nfss,x,ggss;
@#endif
@#if ttb_spec==12
parameters rdscale,tfpwt,tau,  hb,  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, psi, epsilon,   alfa, rhoum,rhouc, rhouh, rhou,rhoud,rhouf, sigc, sigh, sigm, ss, rho_eta, eta1, eta2, eta3, eta4,eta5,eta6,eta7,eta8,eta9,eta10,eta11,eta12,iss, rss, yss, kss, css,wss,mss, rdss, pfss,  hss, ivss, lamdass, nfss,x,ggss;
@#endif
@#if ttb_spec==16
parameters rdscale,tfpwt,tau,  hb,  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, psi, epsilon,   alfa, rhoum,rhouc, rhouh, rhou,rhoud,rhouf, sigc, sigh, sigm, ss, rho_eta, eta1, eta2, eta3, eta4,eta5,eta6,eta7,eta8,eta9,eta10,eta11,eta12,eta13, eta14,eta15,eta16,iss, rss, yss, kss, css,wss,mss, rdss, pfss,  hss, ivss, lamdass, nfss,x,ggss;
@#endif
@#if ttb_spec==20
parameters rdscale,tfpwt,tau,  hb,  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, psi, epsilon,   alfa, rhoum,rhouc, rhouh, rhou,rhoud,rhouf, sigc, sigh, sigm, ss, rho_eta, eta1, eta2, eta3, eta4,eta5,eta6,eta7,eta8,eta9,eta10,eta11,eta12,eta13, eta14,eta15,eta16,eta17,eta18,eta19,eta20,iss, rss, yss, kss, css,wss,mss, rdss, pfss,  hss, ivss, lamdass, nfss,x,ggss;
@#endif
@#if ttb_spec==24
parameters rdscale,tfpwt,tau,  hb,  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, psi, epsilon,   alfa, rhoum,rhouc, rhouh, rhou,rhoud,rhouf, sigc, sigh, sigm, ss, rho_eta, eta1, eta2, eta3, eta4,eta5,eta6,eta7,eta8,eta9,eta10,eta11,eta12,eta13, eta14,eta15,eta16,eta17,eta18,eta19,eta20,eta21,eta22,eta23, eta24,iss, rss, yss, kss, css,wss,mss, rdss, pfss,  hss, ivss, lamdass, nfss,x,ggss;
@#endif
@#if ttb_spec==10
parameters rdscale,tfpwt,tau,  hb,  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, psi, epsilon,   alfa, rhoum,rhouc, rhouh, rhou,rhoud,rhouf, sigc, sigh, sigm, ss, rho_eta, eta1, eta2, eta3, eta4,eta5,eta6,eta7,eta8,eta9,eta10,iss, rss, yss, kss, css,wss,mss, rdss, pfss,  hss, ivss, lamdass, nfss,x,ggss;
@#endif
 
 
@#endif

@#if tech==0
 @#if var_data==20
 var y_g, c_g, iv_g, h_g, w_g, tfp_g, gg, ggn, uf, pai,y,tfp, tfpe, tfpu,u,conf, conf12, conf11, conf10, conf9, conf8,c,k,iv,h,w,r,rr,p,ppai,m,i,lamda,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,in,lamdan,yg,hg,mcgap,uc,um,uh, ynobs, fpai,mcgape, mc,mcn,mce, fpaie, paie, paiee,ppaie,  tfpn,tq,tqn,shu, shi, shc, shh, shue1 sp  stk ;
 @#endif
@#if var_data==21
 var y_g, c_g, iv_g, h_g, w_g, tfp_g, gg, ggn, uf, pai,y,tfp, tfpu,u,conf, conf12, conf11, conf10, conf9, conf8,c,k,iv,h,w,r,rr,p,ppai,m,i,lamda,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,in,lamdan,yg,hg,mcgap,uc,um,uh, ynobs, fpai,mcgape, mc,mcn,mce, fpaie, paie, paiee,ppaie,  tfpn,tq,tqn,shu, shi, shc, shh, shue1 sp  stk ;
 @#endif
@#if var_data==22
 var y_g, c_g, iv_g, h_g, w_g, tfp_g, gg, ggn, uf, pai,y,tfp,tfpe, tfpu,u,conf, conf12, conf11, conf10, conf9, conf8,c,k,iv,h,w,r,rr,p,ppai,m,i,lamda,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,in,lamdan,yg,hg,mcgap,uc,um,uh, ynobs, fpai,mcgape, mc,mcn,mce, fpaie, paie, paiee,ppaie,  tfpn,tq,tqn,shu, shi, shc, shh, shue1 sp  stk ;
 @#endif
@#if var_data==25
 var y_g, c_g, iv_g, h_g, w_g, tfp_g, gg, ggn, uf, pai,y,tfp, tfpu,u,conf, conf12, conf11, conf10, conf9, conf8,c,k,iv,h,w,r,rr,p,ppai,m,i,lamda,ud,v,pain,yn,kn,ivn,hn,cn,wn,rn,pn,in,lamdan,yg,hg,mcgap,uc,um,uh, ynobs, fpai,mcgape, mc,mcn,mce, fpaie, paie, paiee,ppaie,  tfpn,tq,tqn,shu, shi, shc, shh, shue1 sp  stk ;
 @#endif

@#if eshock==16
varexo ei eu euc euh ef eue1 eue2 eue3 eue4 eue5 eue6 eue7 eue8 eue9 eue10 eue11 eue12 eue13 eue14 eue15 eue etfp me mee meee;
@#endif
 @#if pindex==1
 parameters index, tfpwt, tau, hb,  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, rhoum,rhouc,rhouh, rhou,rhoud,rhouf, sigc, sigh, sigm, ss, iss_noa, rss_noa, yss_noa, kss_noa, css_noa,wss_noa,mss_noa,  hss_noa, ivss_noa, lamdass_noa, ggss_noa ;
 @#endif
 @#if pindex==0
 parameters tfpwt, tau, hb,  xxx, rhoi,xpai, xy, gamma, delta, phi, theta, rho, rhoum,rhouc,rhouh, rhou,rhoud,rhouf, sigc, sigh, sigm, ss, iss_noa, rss_noa, yss_noa, kss_noa, css_noa,wss_noa,mss_noa,  hss_noa, ivss_noa, lamdass_noa, ggss_noa ;
 @#endif
@#endif


//** params adjustment**//
rhoaj=rhoaj_;
tfpwt=tfpwt_;
alfaj=alfaj_;
phiaj=phiaj_;
psiaj=psiaj_;
rdscale=rdscale_;
rd_rt=rd_rt_;
//*******************//

@#if tech==1

//** data_techzero==1 section **//

@#if data_techzero==1



@#if eshock==16

 rho_eta=aa0(j,17);
 xxx = aa0(j,18);
 xpai = aa0(j,19); //*1.5 (CLE&Gali) *//
 xy =aa0(j,20); //* 0.5/4 (gali) 0.1 (CLE) *//


 
 @#if calvo_pra==0
 rho =aa0(j,22);  //** 0.75 free param& rho=0 is flexible price
 @#endif

 @#if calvo_pra==1
 rho = rhoaj*aa0(j,22);  //** 0.75 free param& rho=0 is flexible price
 @#endif

 sigc =aa0(j,23); 
 sigm = 2.56; //* CEL=9
 //ss=aa0(j,24);

 @#if simple==0
  hb=aa0(j,21); % Internal habit formation param.
  ss=aa0(j,24);
 @#endif

 @#if simple==2
  hb=0.000001*aa0(j,21); % Internal habit formation param.
  ss=0.000001*aa0(j,24);
 @#endif

 phi = phiaj*aa0(j,25);
 psi = psiaj*aa0(j,26);

 rhou = aa0(j,27);
 rhoi = aa0(j,28);//*0.5 (Gali) &0.15(Smet) *//
 rhoud =0.5;
 rhouh = aa0(j,29); //**this needs to be high to get good results
 rhouc = aa0(j,30); //**this needs to be high to get good results
 rhouf = aa0(j,31); //**this needs to be high to get good results


 @#if ashock==1
 alfa =alfaj*aa0(j,32);
 @#endif
 @#if ashock==0
 alfa =alfa_n;
 @#endif
@#endif


rhoum =0.0;
//sigc=1.00;
gamma =0.99; % discount factor  
delta = 0.10000000/4;  %capital depreciation rate
theta = 1/3;
notec=0;
epsilon = 0.1/4;
//epsilon = 0.00001;
tau =0.18;
//phi=6.0;
sigh = 2.0;

@#if ttb_spec==8
eta1=( (1+rho_eta+(rho_eta^2)+(rho_eta^3)+(rho_eta^4)+(rho_eta^5)+(rho_eta^6)+(rho_eta^7))^(-1) );
eta2=eta1*(rho_eta);
eta3=eta1*(rho_eta)^2; 
eta4=eta1*(rho_eta)^3;
eta5=eta1*(rho_eta)^4;
eta6=eta1*(rho_eta)^5;
eta7=eta1*(rho_eta)^6;
eta8=eta1*(rho_eta)^7;
@#endif

@#if ttb_spec==10
eta1=( (1+rho_eta+(rho_eta^2)+(rho_eta^3)+(rho_eta^4)+(rho_eta^5)+(rho_eta^6)+(rho_eta^7)+(rho_eta^8)+(rho_eta^9))^(-1) );
eta2=eta1*(rho_eta);
eta3=eta1*(rho_eta)^2; 
eta4=eta1*(rho_eta)^3;
eta5=eta1*(rho_eta)^4;
eta6=eta1*(rho_eta)^5;
eta7=eta1*(rho_eta)^6;
eta8=eta1*(rho_eta)^7;
eta9=eta1*(rho_eta)^8;
eta10=eta1*(rho_eta)^9;
@#endif

@#if ttb_spec==12
eta1=( (1+rho_eta+(rho_eta^2)+(rho_eta^3)+(rho_eta^4)+(rho_eta^5)+(rho_eta^6)+(rho_eta^7)+(rho_eta^8)+(rho_eta^9)+(rho_eta^10)+(rho_eta^11))^(-1) );
eta2=eta1*(rho_eta);
eta3=eta1*(rho_eta)^2; 
eta4=eta1*(rho_eta)^3;
eta5=eta1*(rho_eta)^4;
eta6=eta1*(rho_eta)^5;
eta7=eta1*(rho_eta)^6;
eta8=eta1*(rho_eta)^7;
eta9=eta1*(rho_eta)^8;
eta10=eta1*(rho_eta)^9;
eta11=eta1*(rho_eta)^10;
eta12=eta1*(rho_eta)^11;
@#endif


@#if ttb_spec==16
eta1=( (1+rho_eta+(rho_eta^2)+(rho_eta^3)+(rho_eta^4)+(rho_eta^5)+(rho_eta^6)+(rho_eta^7)+(rho_eta^8)+(rho_eta^9)+(rho_eta^10)+(rho_eta^11)+(rho_eta^12)+(rho_eta^13)+(rho_eta^14)+(rho_eta^15))^(-1) );
eta2=eta1*(rho_eta);
eta3=eta1*(rho_eta)^2; 
eta4=eta1*(rho_eta)^3;
eta5=eta1*(rho_eta)^4;
eta6=eta1*(rho_eta)^5;
eta7=eta1*(rho_eta)^6;
eta8=eta1*(rho_eta)^7;
eta9=eta1*(rho_eta)^8;
eta10=eta1*(rho_eta)^9;
eta11=eta1*(rho_eta)^10;
eta12=eta1*(rho_eta)^11;
eta13=eta1*(rho_eta)^12; 
eta14=eta1*(rho_eta)^13;
eta15=eta1*(rho_eta)^14;
eta16=eta1*(rho_eta)^15;
@#endif


@#if ttb_spec==20
eta1=( (1+rho_eta+(rho_eta^2)+(rho_eta^3)+(rho_eta^4)+(rho_eta^5)+(rho_eta^6)+(rho_eta^7)+(rho_eta^8)+(rho_eta^9)+(rho_eta^10)+(rho_eta^11)+(rho_eta^12)+(rho_eta^13)+(rho_eta^14)+(rho_eta^15)+(rho_eta^16)+(rho_eta^17)+(rho_eta^18)+(rho_eta^19)+(rho_eta^20))^(-1) );
eta2=eta1*(rho_eta);
eta3=eta1*(rho_eta)^2; 
eta4=eta1*(rho_eta)^3;
eta5=eta1*(rho_eta)^4;
eta6=eta1*(rho_eta)^5;
eta7=eta1*(rho_eta)^6;
eta8=eta1*(rho_eta)^7;
eta9=eta1*(rho_eta)^8;
eta10=eta1*(rho_eta)^9;
eta11=eta1*(rho_eta)^10;
eta12=eta1*(rho_eta)^11;
eta13=eta1*(rho_eta)^12; 
eta14=eta1*(rho_eta)^13;
eta15=eta1*(rho_eta)^14;
eta16=eta1*(rho_eta)^15;
eta17=eta1*(rho_eta)^16;
eta18=eta1*(rho_eta)^17;
eta19=eta1*(rho_eta)^18;
eta20=eta1*(rho_eta)^19;
@#endif



@#if ttb_spec==24
eta1=( (1+rho_eta+(rho_eta^2)+(rho_eta^3)+(rho_eta^4)+(rho_eta^5)+(rho_eta^6)+(rho_eta^7)+(rho_eta^8)+(rho_eta^9)+(rho_eta^10)+(rho_eta^11)+(rho_eta^12)+(rho_eta^13)+(rho_eta^14)+(rho_eta^15)+(rho_eta^16)+(rho_eta^17)+(rho_eta^18)+(rho_eta^19)+(rho_eta^20)+(rho_eta^21)+(rho_eta^22)+(rho_eta^23)+(rho_eta^24))^(-1) );
eta2=eta1*(rho_eta);
eta3=eta1*(rho_eta)^2; 
eta4=eta1*(rho_eta)^3;
eta5=eta1*(rho_eta)^4;
eta6=eta1*(rho_eta)^5;
eta7=eta1*(rho_eta)^6;
eta8=eta1*(rho_eta)^7;
eta9=eta1*(rho_eta)^8;
eta10=eta1*(rho_eta)^9;
eta11=eta1*(rho_eta)^10;
eta12=eta1*(rho_eta)^11;
eta13=eta1*(rho_eta)^12; 
eta14=eta1*(rho_eta)^13;
eta15=eta1*(rho_eta)^14;
eta16=eta1*(rho_eta)^15;
eta17=eta1*(rho_eta)^16;
eta18=eta1*(rho_eta)^17;
eta19=eta1*(rho_eta)^18;
eta20=eta1*(rho_eta)^19;
eta21=eta1*(rho_eta)^20;
eta22=eta1*(rho_eta)^21;
eta23=eta1*(rho_eta)^22;
eta24=eta1*(rho_eta)^23;
@#endif

iss = (1/gamma)-1;
rss = iss + delta;
wss = (((phi-1)/phi)*(1-theta)*(((1-theta)/theta)^(-theta))*(rss^(-theta)))^(1/(1-theta));
yss_hss = ((theta/(1-theta))*(wss/rss))^theta;
kss_hss = (yss_hss)^(1/theta);
nfss= (1-psi)/epsilon;


@#if ttb_spec==8
g1_t8 = (1/(1-psi*gamma))*(1-psi)*( ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6)  +eta7*((1+iss)^7)  +eta8*((1+iss)^8)  )^(-1) )*(  1-  (1/(1-theta))*( ( (1-theta)/theta )^theta )*( rss^theta)*( wss^(1-theta) )  );  
//g2_t8 = (1-g1_t8)*(yss_hss) - (delta*kss_hss);
g2_t8 = (1-g1_t8)*(yss_hss) - (delta*kss_hss)- tau*(yss_hss);
g1 = g1_t8; 
g2 = g2_t8;
css_hss = g2_t8; 
etas = ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6) +eta7*((1+iss)^7) +eta8*((1+iss)^8) );
@#endif

@#if ttb_spec==10
g1_t10 = (1/(1-psi*gamma))*(1-psi)*( ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6)  +eta7*((1+iss)^7)  +eta8*((1+iss)^8) +eta9*((1+iss)^9) +eta10*((1+iss)^10)     )^(-1) )*(  1-  (1/(1-theta))*( ( (1-theta)/theta )^theta )*( rss^theta)*( wss^(1-theta) )  );  
g2_t10 = (1-g1_t10)*(yss_hss) - (delta*kss_hss);
g1 = g1_t10; 
g2 = g2_t10;
css_hss = g2_t10; 
etas = ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6) +eta7*((1+iss)^7) +eta8*((1+iss)^8)  +eta9*((1+iss)^9)  +eta10*((1+iss)^10)  );
@#endif

@#if ttb_spec==12
g1_t12 = (1/(1-psi*gamma))*(1-psi)*( ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6)  +eta7*((1+iss)^7)  +eta8*((1+iss)^8) +eta9*((1+iss)^9) +eta10*((1+iss)^10)  +eta11*((1+iss)^11)  +eta12*((1+iss)^12)  )^(-1) )*(  1-  (1/(1-theta))*( ( (1-theta)/theta )^theta )*( rss^theta)*( wss^(1-theta) )  );  
g2_t12 = (1-g1_t12)*(yss_hss) - (delta*kss_hss);
g1 = g1_t12; 
g2 = g2_t12;
css_hss = g2_t12; 
etas = ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6) +eta7*((1+iss)^7) +eta8*((1+iss)^8)  +eta9*((1+iss)^9)  +eta10*((1+iss)^10)  +eta11*((1+iss)^11)  +eta12*((1+iss)^12) );
@#endif

@#if ttb_spec==16
g1_t16 = (1/(1-psi*gamma))*(1-psi)*( ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6)  +eta7*((1+iss)^7)  +eta8*((1+iss)^8) +eta9*((1+iss)^9) +eta10*((1+iss)^10)  +eta11*((1+iss)^11)  +eta12*((1+iss)^12) +eta13*((1+iss)^13)+eta14*((1+iss)^14)+eta15*((1+iss)^15)+eta16*((1+iss)^16) )^(-1) )*(  1-  (1/(1-theta))*( ( (1-theta)/theta )^theta )*( rss^theta)*( wss^(1-theta) )  );  
g2_t16 = (1-g1_t16)*(yss_hss) - (delta*kss_hss);
g1 = g1_t16; 
g2 = g2_t16;
css_hss = g2_t16; 
etas = ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6) +eta7*((1+iss)^7) +eta8*((1+iss)^8)  +eta9*((1+iss)^9)  +eta10*((1+iss)^10)  +eta11*((1+iss)^11)  +eta12*((1+iss)^12)+eta13*((1+iss)^13)+eta14*((1+iss)^14)+eta15*((1+iss)^15)+eta16*((1+iss)^16) );
@#endif


@#if ttb_spec==20
g1_t20 = (1/(1-psi*gamma))*(1-psi)*( ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6)  +eta7*((1+iss)^7)  +eta8*((1+iss)^8) +eta9*((1+iss)^9) +eta10*((1+iss)^10)  +eta11*((1+iss)^11)  +eta12*((1+iss)^12) +eta13*((1+iss)^13)+eta14*((1+iss)^14)+eta15*((1+iss)^15)+eta16*((1+iss)^16)+eta17*((1+iss)^17)+eta18*((1+iss)^18)+eta19*((1+iss)^19)+eta20*((1+iss)^20) )^(-1) )*(  1-  (1/(1-theta))*( ( (1-theta)/theta )^theta )*( rss^theta)*( wss^(1-theta) )  );  
g2_t20 = (1-g1_t20)*(yss_hss) - (delta*kss_hss);
g1 = g1_t20; 
g2 = g2_t20;
css_hss = g2_t20; 
etas = ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6) +eta7*((1+iss)^7) +eta8*((1+iss)^8)  +eta9*((1+iss)^9)  +eta10*((1+iss)^10)  +eta11*((1+iss)^11)  +eta12*((1+iss)^12)+eta13*((1+iss)^13)+eta14*((1+iss)^14)+eta15*((1+iss)^15)+eta16*((1+iss)^16) +eta17*((1+iss)^17)+eta18*((1+iss)^18)+eta19*((1+iss)^19)+eta20*((1+iss)^20)  );
@#endif

@#if ttb_spec==24
g1_t24 = (1/(1-psi*gamma))*(1-psi)*( ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6)  +eta7*((1+iss)^7)  +eta8*((1+iss)^8) +eta9*((1+iss)^9) +eta10*((1+iss)^10)  +eta11*((1+iss)^11)  +eta12*((1+iss)^12) +eta13*((1+iss)^13)+eta14*((1+iss)^14)+eta15*((1+iss)^15)+eta16*((1+iss)^16)+eta17*((1+iss)^17)+eta18*((1+iss)^18)+eta19*((1+iss)^19)+eta20*((1+iss)^20) +eta21*((1+iss)^21)+eta22*((1+iss)^22)+eta23*((1+iss)^23)+eta24*((1+iss)^24))^(-1) )*(  1-  (1/(1-theta))*( ( (1-theta)/theta )^theta )*( rss^theta)*( wss^(1-theta) )  );  
g2_t24 = (1-g1_t24)*(yss_hss) - (delta*kss_hss);
g1 = g1_t24; 
g2 = g2_t24;
css_hss = g2_t24; 
etas = ( eta1*(1+iss) + eta2*((1+iss)^2) + eta3*((1+iss)^3) +eta4*((1+iss)^4) +eta5*((1+iss)^5) +eta6*((1+iss)^6) +eta7*((1+iss)^7) +eta8*((1+iss)^8)  +eta9*((1+iss)^9)  +eta10*((1+iss)^10)  +eta11*((1+iss)^11)  +eta12*((1+iss)^12)+eta13*((1+iss)^13)+eta14*((1+iss)^14)+eta15*((1+iss)^15)+eta16*((1+iss)^16) +eta17*((1+iss)^17)+eta18*((1+iss)^18)+eta19*((1+iss)^19)+eta20*((1+iss)^20) +eta21*((1+iss)^21)+eta22*((1+iss)^22)+eta23*((1+iss)^23)+eta24*((1+iss)^24) );
@#endif

hss= (wss*(css_hss^(-sigc))*((1-hb)^(-sigc))*(1-hb*gamma))^(1/(sigh+sigc)); 
css=(css_hss)*hss;
yss = (((theta/(1-theta))*(wss/rss))^theta)*hss;
kss = (yss_hss^(1/theta))*hss;

rdss =rd_rt*g1*yss;
//rdss = g1*yss;
x=((nfss)^(1+alfa))/rdss; % x is (kappa_bar*delta) in the paper
//x=((nfss)^(1+0.2509))/rdss; % x is (kappa_bar*delta) in the paper

pfss = (1/(1-psi))*etas*rdss; % pfss is the steady state value of th epresent value of monopoly profits stream 
lamdass = (1-hb*gamma)*((1-hb)*css)^(-sigc); % lamdass = the steady state velu of CHI in the paper
mss = (((1-hb)*css)^(sigc/sigm))*(((1+iss)/iss)^(1/sigm))*((1-hb*gamma)^(-1/sigm));
ivss = delta*kss;
ggss =tau*yss; % ggss=staedy state value of gg, govt expnditure

@#endif
//** End of data_techzero==1 section **//


@#endif




@#if tech==0

@#if data_techzero==0  

@#if eshock==16
tfpwt=tfpwt_;

 @#if pindex==1
 index=aa0(j,17); 
 xxx = aa0(j,18);
 xpai = aa0(j,19); //*1.5 (CLE&Gali) *//
 xy =aa0(j,20); //* 0.5/4 (gali) 0.1 (CLE) *//

//hb=aa0(j,21); % Internal habit formation param.
rho =aa0(j,22);  //** 0.75 free param& rho=0 is flexible price
//sigh =aa0(j,23); 
sigc =aa0(j,23); 
sigm = 2.56; //* CEL=9
//ss=aa0(j,24);
phi = phiaj*aa0(j,25);

@#if simple==0
  hb=aa0(j,21); % Internal habit formation param.
  ss=aa0(j,24);
@#endif

@#if simple==2
  hb=0.000001*aa0(j,21); % Internal habit formation param.
  ss=0.000001*aa0(j,24);
@#endif


 rhou = aa0(j,26);
 rhoi = aa0(j,27);//*0.5 (Gali) &0.15(Smet) *//
 rhoud = 0.5;
 rhouh = aa0(j,28); //**this needs to be high to get good results
 rhouc = aa0(j,29); //**this needs to be high to get good results
 rhouf = aa0(j,30); //**this needs to be high to get good results
 @#endif 
 
 @#if pindex==0
xxx = aa0(j,17);
xpai = aa0(j,18); //*1.5 (CLE&Gali) *//
xy =aa0(j,19); //* 0.5/4 (gali) 0.1 (CLE) *//

//hb=aa0(j,20); % Internal habit formation param.
rho =aa0(j,21);  //** 0.75 free param& rho=0 is flexible price
//sigh =aa0(j,22); 
sigc =aa0(j,22); 
sigm = 2.56; //* CEL=9
//ss=aa0(j,23);
phi = phiaj*aa0(j,24);

@#if simple==0
  hb=aa0(j,20); % Internal habit formation param.
  ss=aa0(j,23);
@#endif

@#if simple==2
  hb=0.000001*aa0(j,20); % Internal habit formation param.
  ss=0.000001*aa0(j,23);
@#endif

rhou = aa0(j,25);
rhoi = aa0(j,26);//*0.5 (Gali) &0.15(Smet) *//
rhoud = 0.5;
rhouh = aa0(j,27); //**this needs to be high to get good results
rhouc = aa0(j,28); //**this needs to be high to get good results
rhouf = aa0(j,29); //**this needs to be high to get good results
@#endif

@#endif


rhoum =0.9;
//sigc=1;
gamma =0.99; % discount factor  
delta = 0.10000000/4;  %capital depreciation rate
theta = 1/3;
notec=0;
tau=0.18;
//phi=6.0;
sigh=2;

iss_noa = (1/gamma)-1;
rss_noa = iss_noa + delta;
wss_noa = (((phi-1)/phi)*(1-theta)*(((1-theta)/theta)^(-theta))*(rss_noa^(-theta)))^(1/(1-theta));
yss_hss_noa = ((theta/(1-theta))*(wss_noa/rss_noa))^theta;
kss_hss_noa = (yss_hss_noa)^(1/theta);
css_hss_noa = yss_hss_noa - delta*kss_hss_noa- tau*yss_hss_noa; //* including the gov expenditure part "-tau*yss_hss_noa"
hss_noa = (wss_noa*(css_hss_noa^(-sigc))*((1-hb)^(-sigc))*(1-hb*gamma))^(1/(sigh+sigc));
yss_noa = (((theta/(1-theta))*(wss_noa/rss_noa))^theta)*hss_noa;
kss_noa = (yss_hss_noa^(1/theta))*hss_noa;
css_noa = css_hss_noa*hss_noa;
mss_noa = (((1-hb)*css_noa)^(sigc/sigm))*(((1+iss_noa)/iss_noa)^(1/sigm))*((1-hb*gamma)^(-1/sigm));
ivss_noa = delta*kss_noa;
ggss_noa =tau*yss_noa; % ggss=staedy state value of gg, govt expnditure

lamdass_noa = (1-hb*gamma)*((1-hb)*css_noa)^(-sigc);
@#endif

@#endif




model(linear);

@#if tech==1


@#if ttb_spec==8

 //rd=(1+alfa)*(eta1*nf(-7)+eta2*nf(-6)+eta3*nf(-5)+eta4*nf(-4)+eta5*nf(-3)+eta6*nf(-2)+eta7*nf(-1)+eta8*nf);
 //pf(+8)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8)*(1/iss)*i(+7) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8)*(1/iss)*i(+6) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8)*(1/iss)*i(+5) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8)*(1/iss)*i(+4) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8)*(1/iss)*i(+3) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8)*(1/iss)*i(+2) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8)*(1/iss)*i(+1) + iss*((1+iss)^7)*(eta8)*(1/iss)*i + eta1*(1+iss)*(p(+7)-p(+8)) + eta2*((1+iss)^2)*(p(+6)-p(+8)) + eta3*((1+iss)^3)*(p(+5)-p(+8)) + eta4*((1+iss)^4)*(p(+4)-p(+8)) + eta5*((1+iss)^5)*(p(+3) -p(+8)) + eta6*((1+iss)^6)*(p(+2) -p(+8)) + eta7*((1+iss)^7)*(p(+1) -p(+8)) + eta8*((1+iss)^8)*(p -p(+8))  );    
 //a = (1-psi)*nf(-7) + psi*a(-1);

 @#if TI==1
 rd=(1+alfa)*(eta1*nf(-7)+eta2*nf(-6)+eta3*nf(-5)+eta4*nf(-4)+eta5*nf(-3)+eta6*nf(-2)+eta7*nf(-1)+eta8*nf);
 pf(+8)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8)*(1/iss)*i(+7) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8)*(1/iss)*i(+6) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8)*(1/iss)*i(+5) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8)*(1/iss)*i(+4) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8)*(1/iss)*i(+3) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8)*(1/iss)*i(+2) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8)*(1/iss)*i(+1) + iss*((1+iss)^7)*(eta8)*(1/iss)*i 
 + eta1*(1+iss)*(p(+7)-p(+8)) + eta2*((1+iss)^2)*(p(+6)-p(+8)) + eta3*((1+iss)^3)*(p(+5)-p(+8)) + eta4*((1+iss)^4)*(p(+4)-p(+8)) + eta5*((1+iss)^5)*(p(+3) -p(+8)) + eta6*((1+iss)^6)*(p(+2) -p(+8)) + eta7*((1+iss)^7)*(p(+1) -p(+8)) + eta8*((1+iss)^8)*(p -p(+8))  );    
 a = rdscale*(1-psi)*nf(-7) + psi*a(-1);
 @#endif

  @#if TI==0
 rd=(1+alfa)*(nf);
 pf(+1)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* ( (i)  +(1+iss)*(p-p(+1)) );    
 a = rdscale*(1-psi)*nf + psi*a(-1);
 @#endif

@#endif

 
@#if ttb_spec==12
 @#if TI==12

 rd=(1+alfa)*(eta1*nf(-11)+eta2*nf(-10)+eta3*nf(-9)+eta4*nf(-8)+eta5*nf(-7)+eta6*nf(-6)+eta7*nf(-5)+eta8*nf(-4)+eta9*nf(-3)+eta10*nf(-2)+eta11*nf(-1)+eta12*nf);
 pf(+12)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10 +((1+iss)^10)*eta11 +((1+iss)^11)*eta12)*(1/iss)*i(+11) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10 +((1+iss)^9)*eta11 +((1+iss)^10)*eta12)*(1/iss)*i(+10) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10 +((1+iss)^8)*eta11 +((1+iss)^9)*eta12)*(1/iss)*i(+9) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10 +((1+iss)^7)*eta11 +((1+iss)^8)*eta12)*(1/iss)*i(+8) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10 +((1+iss)^6)*eta11 +((1+iss)^7)*eta12)*(1/iss)*i(+7) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10 +((1+iss)^5)*eta11 +((1+iss)^6)*eta12)*(1/iss)*i(+6) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10 +((1+iss)^4)*eta11+((1+iss)^5)*eta12)*(1/iss)*i(+5) + iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10 +((1+iss)^3)*eta11 +((1+iss)^4)*eta12)*(1/iss)*i(+4) + iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10 +((1+iss)^2)*eta11 +((1+iss)^3)*eta12)*(1/iss)*i(+3) + iss*((1+iss)^9)*(eta10 +((1+iss)^1)*eta11 +((1+iss)^2)*eta12)*(1/iss)*i(+2)  + iss*((1+iss)^10)*(eta11 +((1+iss)^1)*eta12)*(1/iss)*i(+1)  + iss*((1+iss)^11)*(eta12)*(1/iss)*i + eta1*(1+iss)*(p(+11)-p(+12)) + eta2*((1+iss)^2)*(p(+10)-p(+12)) + eta3*((1+iss)^3)*(p(+9)-p(+12)) + eta4*((1+iss)^4)*(p(+8)-p(+12)) + eta5*((1+iss)^5)*(p(+7) -p(+12)) + eta6*((1+iss)^6)*(p(+6) -p(+12)) + eta7*((1+iss)^7)*(p(+5) -p(+12)) + eta8*((1+iss)^8)*(p(+4) -p(+12)) + eta9*((1+iss)^9)*(p(+3) -p(+12)) + eta10*((1+iss)^10)*(p(+2) -p(+12)) + eta11*((1+iss)^11)*(p(+1) -p(+12)) + eta12*((1+iss)^12)*(p -p(+12))  );    
 a = rdscale*(1-psi)*nf(-11) + psi*a(-1);
  //rd=(1+alfa)*(eta1*nf(-15)+eta2*nf(-14)+eta3*nf(-13)+eta4*nf(-12)+eta5*nf(-11)+eta6*nf(-10)+eta7*nf(-9)+eta8*nf);
 //pf(+16)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8)*(1/iss)*i(+7) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8)*(1/iss)*i(+6) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8)*(1/iss)*i(+5) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8)*(1/iss)*i(+4) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8)*(1/iss)*i(+3) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8)*(1/iss)*i(+2) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8)*(1/iss)*i(+1) + iss*((1+iss)^7)*(eta8)*(1/iss)*i + eta1*(1+iss)*(p(+7)-p(+8)) + eta2*((1+iss)^2)*(p(+6)-p(+8)) + eta3*((1+iss)^3)*(p(+5)-p(+8)) + eta4*((1+iss)^4)*(p(+4)-p(+8)) + eta5*((1+iss)^5)*(p(+3) -p(+8)) + eta6*((1+iss)^6)*(p(+2) -p(+8)) + eta7*((1+iss)^7)*(p(+1) -p(+8)) + eta8*((1+iss)^8)*(p -p(+8))  );    
 //a = (1-psi)*nf(-7) + psi*a(-1);
 @#endif
@#endif


@#if ttb_spec==16
 @#if TI==16

 rd=(1+alfa)*(eta1*nf(-15)+eta2*nf(-14)+eta3*nf(-13)+eta4*nf(-12)+eta5*nf(-11)+eta6*nf(-10)+eta7*nf(-9)+eta8*nf(-8)+eta9*nf(-7)+eta10*nf(-6)+eta11*nf(-5)+eta12*nf*(-4)+eta13*nf*(-3)+eta14*nf*(-2)+eta15*nf*(-1)+eta16*nf);
 pf(+16)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* 
(  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10 +((1+iss)^10)*eta11 +((1+iss)^11)*eta12+((1+iss)^12)*eta13+((1+iss)^13)*eta14+((1+iss)^14)*eta15+((1+iss)^15)*eta16)*(1/iss)*i(+15) 
+  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10 +((1+iss)^9)*eta11 +((1+iss)^10)*eta12 +((1+iss)^11)*eta13 +((1+iss)^12)*eta14 +((1+iss)^13)*eta15 +((1+iss)^14)*eta16)*(1/iss)*i(+14) 
+  iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10 +((1+iss)^8)*eta11 +((1+iss)^9)*eta12 +((1+iss)^10)*eta13  +((1+iss)^11)*eta14 +((1+iss)^12)*eta15 +((1+iss)^13)*eta16)*(1/iss)*i(+13) 
+  iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10 +((1+iss)^7)*eta11 +((1+iss)^8)*eta12 +((1+iss)^9)*eta13 +((1+iss)^10)*eta14 +((1+iss)^11)*eta15 +((1+iss)^12)*eta16 )*(1/iss)*i(+12) 
+  iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10 +((1+iss)^6)*eta11 +((1+iss)^7)*eta12 +((1+iss)^8)*eta13 +((1+iss)^9)*eta14 +((1+iss)^10)*eta15 +((1+iss)^11)*eta16)*(1/iss)*i(+11) 
+  iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10 +((1+iss)^5)*eta11 +((1+iss)^6)*eta12 +((1+iss)^7)*eta13 +((1+iss)^8)*eta14 +((1+iss)^9)*eta15 +((1+iss)^10)*eta16)*(1/iss)*i(+10) 
+  iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10 +((1+iss)^4)*eta11+((1+iss)^5)*eta12 +((1+iss)^6)*eta13 +((1+iss)^7)*eta14 +((1+iss)^8)*eta15 +((1+iss)^9)*eta16)*(1/iss)*i(+9) 
+  iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10 +((1+iss)^3)*eta11 +((1+iss)^4)*eta12  +((1+iss)^5)*eta13  +((1+iss)^6)*eta14  +((1+iss)^7)*eta15  +((1+iss)^8)*eta16)*(1/iss)*i(+8) 
+  iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10 +((1+iss)^2)*eta11 +((1+iss)^3)*eta12 +((1+iss)^4)*eta13  +((1+iss)^5)*eta14  +((1+iss)^6)*eta15  +((1+iss)^7)*eta16 )*(1/iss)*i(+7) 
+  iss*((1+iss)^9)*(eta10 +((1+iss)^1)*eta11 +((1+iss)^2)*eta12 +((1+iss)^3)*eta13 +((1+iss)^4)*eta14  +((1+iss)^5)*eta15  +((1+iss)^6)*eta16  )*(1/iss)*i(+6) 
+  iss*((1+iss)^10)*(eta11 +((1+iss)^1)*eta12 +((1+iss)^2)*eta13 +((1+iss)^3)*eta14 +((1+iss)^4)*eta15  +((1+iss)^5)*eta16    )*(1/iss)*i(+5) 
+  iss*((1+iss)^11)*(eta12 +((1+iss)^1)*eta13 +((1+iss)^2)*eta14 +((1+iss)^3)*eta15 +((1+iss)^4)*eta16      )*(1/iss)*i(+4) 
+  iss*((1+iss)^12)*(eta13 +((1+iss)^1)*eta14 +((1+iss)^2)*eta15 +((1+iss)^3)*eta16    )*(1/iss)*i(+3) 
+  iss*((1+iss)^13)*(eta14 +((1+iss)^1)*eta15 +((1+iss)^2)*eta16    )*(1/iss)*i(+2)  
+  iss*((1+iss)^14)*(eta15 +((1+iss)^1)*eta16    )*(1/iss)*i(+1)  
+  iss*((1+iss)^15)*(eta16   )*(1/iss)*i 
+  eta1*(1+iss)*(p(+15)-p(+16)) + eta2*((1+iss)^2)*(p(+14)-p(+16)) + eta3*((1+iss)^3)*(p(+13)-p(+16)) + eta4*((1+iss)^4)*(p(+12)-p(+16)) + eta5*((1+iss)^5)*(p(+11) -p(+16)) + eta6*((1+iss)^6)*(p(+10) -p(+16)) + eta7*((1+iss)^7)*(p(+9) -p(+16)) + eta8*((1+iss)^8)*(p(+8) -p(+16)) + eta9*((1+iss)^9)*(p(+7) -p(+16)) + eta10*((1+iss)^10)*(p(+6) -p(+16)) + eta11*((1+iss)^11)*(p(+5) -p(+16)) + eta12*((1+iss)^12)*(p(+4) -p(+16)) + eta13*((1+iss)^13)*(p(+3) -p(+16))+ eta14*((1+iss)^14)*(p(+2) -p(+16))+ eta15*((1+iss)^15)*(p(+1) -p(+16))+ eta16*((1+iss)^16)*(p -p(+16))  );    
a = rdscale*(1-psi)*nf(-15) + psi*a(-1);
  
 @#endif
@#endif



@#if ttb_spec==20
 @#if TI==20

 rd=(1+alfa)*(eta1*nf(-19)+eta2*nf(-18)+eta3*nf(-17)+eta4*nf(-16)+eta5*nf(-15)+eta6*nf(-14)+eta7*nf(-13)+eta8*nf(-12)+eta9*nf(-11)+eta10*nf(-10)+eta11*nf(-9)+eta12*nf(-8)+eta13*nf(-7)+eta14*nf(-6)+eta15*nf(-5)+eta16*nf*(-4)+eta17*nf*(-3)+eta18*nf*(-2)+eta19*nf*(-1)+eta20*nf);          
 pf(+20)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* 
(  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10 +((1+iss)^10)*eta11 +((1+iss)^11)*eta12+((1+iss)^12)*eta13+((1+iss)^13)*eta14 +((1+iss)^14)*eta15+((1+iss)^15)*eta16+((1+iss)^16)*eta17+((1+iss)^17)*eta18+((1+iss)^18)*eta19 +((1+iss)^19)*eta20  )*(1/iss)*i(+19) 
+  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10 +((1+iss)^9)*eta11 +((1+iss)^10)*eta12 +((1+iss)^11)*eta13 +((1+iss)^12)*eta14 +((1+iss)^13)*eta15 +((1+iss)^14)*eta16  +((1+iss)^15)*eta17 +((1+iss)^16)*eta18 +((1+iss)^17)*eta19 +((1+iss)^18)*eta20  )*(1/iss)*i(+18) 
+  iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10 +((1+iss)^8)*eta11 +((1+iss)^9)*eta12 +((1+iss)^10)*eta13  +((1+iss)^11)*eta14 +((1+iss)^12)*eta15 +((1+iss)^13)*eta16 +((1+iss)^14)*eta17 +((1+iss)^15)*eta18 +((1+iss)^16)*eta19 +((1+iss)^17)*eta20  )*(1/iss)*i(+17) 
+  iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10 +((1+iss)^7)*eta11 +((1+iss)^8)*eta12 +((1+iss)^9)*eta13 +((1+iss)^10)*eta14 +((1+iss)^11)*eta15 +((1+iss)^12)*eta16 +((1+iss)^13)*eta17 +((1+iss)^14)*eta18 +((1+iss)^15)*eta19 +((1+iss)^16)*eta20  )*(1/iss)*i(+16) 
+  iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10 +((1+iss)^6)*eta11 +((1+iss)^7)*eta12 +((1+iss)^8)*eta13 +((1+iss)^9)*eta14 +((1+iss)^10)*eta15 +((1+iss)^11)*eta16 +((1+iss)^12)*eta17 +((1+iss)^13)*eta18 +((1+iss)^14)*eta19 +((1+iss)^15)*eta20  )*(1/iss)*i(+15) 

+  iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10 +((1+iss)^5)*eta11 +((1+iss)^6)*eta12 +((1+iss)^7)*eta13 +((1+iss)^8)*eta14 +((1+iss)^9)*eta15 +((1+iss)^10)*eta16 +((1+iss)^11)*eta17 +((1+iss)^12)*eta18 +((1+iss)^13)*eta19 +((1+iss)^14)*eta20  )*(1/iss)*i(+14) 
+  iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10 +((1+iss)^4)*eta11+((1+iss)^5)*eta12 +((1+iss)^6)*eta13 +((1+iss)^7)*eta14 +((1+iss)^8)*eta15   +((1+iss)^9)*eta16 +((1+iss)^10)*eta17 +((1+iss)^11)*eta18 +((1+iss)^12)*eta19 +((1+iss)^13)*eta20  )*(1/iss)*i(+13) 
+  iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10 +((1+iss)^3)*eta11 +((1+iss)^4)*eta12  +((1+iss)^5)*eta13  +((1+iss)^6)*eta14  +((1+iss)^7)*eta15 +((1+iss)^8)*eta16 +((1+iss)^9)*eta17 +((1+iss)^10)*eta18 +((1+iss)^11)*eta19 +((1+iss)^12)*eta20  )*(1/iss)*i(+12) 
+  iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10 +((1+iss)^2)*eta11 +((1+iss)^3)*eta12 +((1+iss)^4)*eta13  +((1+iss)^5)*eta14  +((1+iss)^6)*eta15  +((1+iss)^7)*eta16 +((1+iss)^8)*eta17 +((1+iss)^9)*eta18 +((1+iss)^10)*eta19 +((1+iss)^11)*eta20   )*(1/iss)*i(+11) 
+  iss*((1+iss)^9)*(eta10 +((1+iss)^1)*eta11 +((1+iss)^2)*eta12 +((1+iss)^3)*eta13 +((1+iss)^4)*eta14  +((1+iss)^5)*eta15  +((1+iss)^6)*eta16 +((1+iss)^7)*eta17 +((1+iss)^8)*eta18 +((1+iss)^9)*eta19 +((1+iss)^10)*eta20    )*(1/iss)*i(+10) 
+  iss*((1+iss)^10)*(eta11 +((1+iss)^1)*eta12 +((1+iss)^2)*eta13 +((1+iss)^3)*eta14 +((1+iss)^4)*eta15  +((1+iss)^5)*eta16 +((1+iss)^6)*eta17 +((1+iss)^7)*eta18 +((1+iss)^8)*eta19 +((1+iss)^9)*eta20     )*(1/iss)*i(+9) 
+  iss*((1+iss)^11)*(eta12 +((1+iss)^1)*eta13 +((1+iss)^2)*eta14 +((1+iss)^3)*eta15 +((1+iss)^4)*eta16 +((1+iss)^5)*eta17 +((1+iss)^6)*eta18 +((1+iss)^7)*eta19 +((1+iss)^8)*eta20     )*(1/iss)*i(+8) 
+  iss*((1+iss)^12)*(eta13 +((1+iss)^1)*eta14 +((1+iss)^2)*eta15 +((1+iss)^3)*eta16 +((1+iss)^4)*eta17 +((1+iss)^5)*eta18 +((1+iss)^6)*eta19 +((1+iss)^7)*eta20    )*(1/iss)*i(+7) 
+  iss*((1+iss)^13)*(eta14 +((1+iss)^1)*eta15 +((1+iss)^2)*eta16 +((1+iss)^3)*eta17 +((1+iss)^4)*eta18 +((1+iss)^5)*eta19 +((1+iss)^6)*eta20    )*(1/iss)*i(+6)  
+  iss*((1+iss)^14)*(eta15 +((1+iss)^1)*eta16 +((1+iss)^2)*eta17 +((1+iss)^3)*eta18 +((1+iss)^4)*eta19 +((1+iss)^5)*eta20   )*(1/iss)*i(+5)  
+  iss*((1+iss)^15)*(eta16 +((1+iss)^1)*eta17 +((1+iss)^2)*eta18 +((1+iss)^3)*eta19 +((1+iss)^4)*eta20  )*(1/iss)*i(+4)  
+  iss*((1+iss)^16)*(eta17 +((1+iss)^1)*eta18 +((1+iss)^2)*eta19 +((1+iss)^3)*eta20    )*(1/iss)*i(+3)  
+  iss*((1+iss)^17)*(eta18 +((1+iss)^1)*eta19 +((1+iss)^2)*eta20   )*(1/iss)*i(+2)  
+  iss*((1+iss)^18)*(eta19 +((1+iss)^1)*eta20   )*(1/iss)*i(+1)  
+  iss*((1+iss)^19)*(eta20     )*(1/iss)*i 
+  eta1*(1+iss)*(p(+19)-p(+20)) + eta2*((1+iss)^2)*(p(+18)-p(+20)) + eta3*((1+iss)^3)*(p(+17)-p(+20)) + eta4*((1+iss)^4)*(p(+16)-p(+20)) + eta5*((1+iss)^5)*(p(+15) -p(+20)) + eta6*((1+iss)^6)*(p(+14) -p(+20)) + eta7*((1+iss)^7)*(p(+13) -p(+20)) + eta8*((1+iss)^8)*(p(+12) -p(+20)) + eta9*((1+iss)^9)*(p(+11) -p(+20)) + eta10*((1+iss)^10)*(p(+10) -p(+20)) + eta11*((1+iss)^11)*(p(+9) -p(+20)) + eta12*((1+iss)^12)*(p(+8) -p(+20)) + eta13*((1+iss)^13)*(p(+7)-p(+20))+ eta14*((1+iss)^14)*(p(+6) -p(+20)) + eta15*((1+iss)^15)*(p(+5) -p(+20))  + eta16*((1+iss)^16)*(p(+4) -p(+20)) + eta17*((1+iss)^17)*(p(+3) -p(+20)) + eta18*((1+iss)^18)*(p(+2) -p(+20)) + eta19*((1+iss)^19)*(p(+1) -p(+20)) + eta20*((1+iss)^20)*(p -p(+20))  );    

a = rdscale*(1-psi)*nf(-19) + psi*a(-1);
  @#endif
@#endif




@#if ttb_spec==24
 @#if TI==24

 rd=(1+alfa)*(eta1*nf(-23)+eta2*nf(-22)+eta3*nf(-21)+eta4*nf(-20)+eta5*nf(-19)+eta6*nf(-18)+eta7*nf(-17)+eta8*nf(-16)+eta9*nf(-15)+eta10*nf(-14)+eta11*nf(-13)+eta12*nf(-12)+eta13*nf(-11)+eta14*nf(-10)+eta15*nf(-9)+eta16*nf(-8)+eta17*nf(-7)+eta18*nf(-6)+eta19*nf(-5)+eta20*nf*(-4)+eta21*nf*(-3)+eta22*nf*(-2)+eta23*nf*(-1)+eta24*nf);          
 pf(+24)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* 
(  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10 +((1+iss)^10)*eta11 +((1+iss)^11)*eta12+((1+iss)^12)*eta13+((1+iss)^13)*eta14 +((1+iss)^14)*eta15+((1+iss)^15)*eta16+((1+iss)^16)*eta17+((1+iss)^17)*eta18+((1+iss)^18)*eta19 +((1+iss)^19)*eta20  +((1+iss)^20)*eta21 +((1+iss)^21)*eta22 +((1+iss)^22)*eta23 +((1+iss)^23)*eta24)*(1/iss)*i(+23) 
+  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10 +((1+iss)^9)*eta11 +((1+iss)^10)*eta12 +((1+iss)^11)*eta13 +((1+iss)^12)*eta14 +((1+iss)^13)*eta15 +((1+iss)^14)*eta16  +((1+iss)^15)*eta17 +((1+iss)^16)*eta18 +((1+iss)^17)*eta19 +((1+iss)^18)*eta20  +((1+iss)^19)*eta21 +((1+iss)^20)*eta22 +((1+iss)^21)*eta23 +((1+iss)^22)*eta24)*(1/iss)*i(+22) 
+  iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10 +((1+iss)^8)*eta11 +((1+iss)^9)*eta12 +((1+iss)^10)*eta13  +((1+iss)^11)*eta14 +((1+iss)^12)*eta15 +((1+iss)^13)*eta16 +((1+iss)^14)*eta17 +((1+iss)^15)*eta18 +((1+iss)^16)*eta19 +((1+iss)^17)*eta20  +((1+iss)^18)*eta21 +((1+iss)^19)*eta22 +((1+iss)^20)*eta23 +((1+iss)^21)*eta24)*(1/iss)*i(+21) 
+  iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10 +((1+iss)^7)*eta11 +((1+iss)^8)*eta12 +((1+iss)^9)*eta13 +((1+iss)^10)*eta14 +((1+iss)^11)*eta15 +((1+iss)^12)*eta16 +((1+iss)^13)*eta17 +((1+iss)^14)*eta18 +((1+iss)^15)*eta19 +((1+iss)^16)*eta20  +((1+iss)^17)*eta21 +((1+iss)^18)*eta22 +((1+iss)^19)*eta23 +((1+iss)^20)*eta24)*(1/iss)*i(+20) 
+  iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10 +((1+iss)^6)*eta11 +((1+iss)^7)*eta12 +((1+iss)^8)*eta13 +((1+iss)^9)*eta14 +((1+iss)^10)*eta15 +((1+iss)^11)*eta16 +((1+iss)^12)*eta17 +((1+iss)^13)*eta18 +((1+iss)^14)*eta19 +((1+iss)^15)*eta20  +((1+iss)^16)*eta21 +((1+iss)^17)*eta22 +((1+iss)^18)*eta23 +((1+iss)^19)*eta24)*(1/iss)*i(+19) 

+  iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10 +((1+iss)^5)*eta11 +((1+iss)^6)*eta12 +((1+iss)^7)*eta13 +((1+iss)^8)*eta14 +((1+iss)^9)*eta15 +((1+iss)^10)*eta16 +((1+iss)^11)*eta17 +((1+iss)^12)*eta18 +((1+iss)^13)*eta19 +((1+iss)^14)*eta20  +((1+iss)^15)*eta21 +((1+iss)^16)*eta22 +((1+iss)^17)*eta23 +((1+iss)^18)*eta24)*(1/iss)*i(+18) 
+  iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10 +((1+iss)^4)*eta11+((1+iss)^5)*eta12 +((1+iss)^6)*eta13 +((1+iss)^7)*eta14 +((1+iss)^8)*eta15   +((1+iss)^9)*eta16 +((1+iss)^10)*eta17 +((1+iss)^11)*eta18 +((1+iss)^12)*eta19 +((1+iss)^13)*eta20  +((1+iss)^14)*eta21 +((1+iss)^15)*eta22 +((1+iss)^16)*eta23 +((1+iss)^17)*eta24)*(1/iss)*i(+17) 
+  iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10 +((1+iss)^3)*eta11 +((1+iss)^4)*eta12  +((1+iss)^5)*eta13  +((1+iss)^6)*eta14  +((1+iss)^7)*eta15 +((1+iss)^8)*eta16 +((1+iss)^9)*eta17 +((1+iss)^10)*eta18 +((1+iss)^11)*eta19 +((1+iss)^12)*eta20  +((1+iss)^13)*eta21 +((1+iss)^14)*eta22 +((1+iss)^15)*eta23 +((1+iss)^16)*eta24)*(1/iss)*i(+16) 
+  iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10 +((1+iss)^2)*eta11 +((1+iss)^3)*eta12 +((1+iss)^4)*eta13  +((1+iss)^5)*eta14  +((1+iss)^6)*eta15  +((1+iss)^7)*eta16 +((1+iss)^8)*eta17 +((1+iss)^9)*eta18 +((1+iss)^10)*eta19 +((1+iss)^11)*eta20  +((1+iss)^12)*eta21 +((1+iss)^13)*eta22 +((1+iss)^14)*eta23 +((1+iss)^15)*eta24 )*(1/iss)*i(+15) 
+  iss*((1+iss)^9)*(eta10 +((1+iss)^1)*eta11 +((1+iss)^2)*eta12 +((1+iss)^3)*eta13 +((1+iss)^4)*eta14  +((1+iss)^5)*eta15  +((1+iss)^6)*eta16 +((1+iss)^7)*eta17 +((1+iss)^8)*eta18 +((1+iss)^9)*eta19 +((1+iss)^10)*eta20  +((1+iss)^11)*eta21 +((1+iss)^12)*eta22 +((1+iss)^13)*eta23 +((1+iss)^14)*eta24  )*(1/iss)*i(+14) 
+  iss*((1+iss)^10)*(eta11 +((1+iss)^1)*eta12 +((1+iss)^2)*eta13 +((1+iss)^3)*eta14 +((1+iss)^4)*eta15  +((1+iss)^5)*eta16 +((1+iss)^6)*eta17 +((1+iss)^7)*eta18 +((1+iss)^8)*eta19 +((1+iss)^9)*eta20  +((1+iss)^10)*eta21 +((1+iss)^11)*eta22 +((1+iss)^12)*eta23 +((1+iss)^13)*eta24    )*(1/iss)*i(+13) 
+  iss*((1+iss)^11)*(eta12 +((1+iss)^1)*eta13 +((1+iss)^2)*eta14 +((1+iss)^3)*eta15 +((1+iss)^4)*eta16 +((1+iss)^5)*eta17 +((1+iss)^6)*eta18 +((1+iss)^7)*eta19 +((1+iss)^8)*eta20  +((1+iss)^9)*eta21 +((1+iss)^10)*eta22 +((1+iss)^11)*eta23 +((1+iss)^12)*eta24      )*(1/iss)*i(+12) 
+  iss*((1+iss)^12)*(eta13 +((1+iss)^1)*eta14 +((1+iss)^2)*eta15 +((1+iss)^3)*eta16 +((1+iss)^4)*eta17 +((1+iss)^5)*eta18 +((1+iss)^6)*eta19 +((1+iss)^7)*eta20  +((1+iss)^8)*eta21 +((1+iss)^9)*eta22 +((1+iss)^10)*eta23 +((1+iss)^11)*eta24    )*(1/iss)*i(+11) 
+  iss*((1+iss)^13)*(eta14 +((1+iss)^1)*eta15 +((1+iss)^2)*eta16 +((1+iss)^3)*eta17 +((1+iss)^4)*eta18 +((1+iss)^5)*eta19 +((1+iss)^6)*eta20  +((1+iss)^7)*eta21 +((1+iss)^8)*eta22 +((1+iss)^9)*eta23 +((1+iss)^10)*eta24    )*(1/iss)*i(+10)  
+  iss*((1+iss)^14)*(eta15 +((1+iss)^1)*eta16 +((1+iss)^2)*eta17 +((1+iss)^3)*eta18 +((1+iss)^4)*eta19 +((1+iss)^5)*eta20  +((1+iss)^6)*eta21 +((1+iss)^7)*eta22 +((1+iss)^8)*eta23 +((1+iss)^9)*eta24    )*(1/iss)*i(+9)  
+  iss*((1+iss)^15)*(eta16 +((1+iss)^1)*eta17 +((1+iss)^2)*eta18 +((1+iss)^3)*eta19 +((1+iss)^4)*eta20  +((1+iss)^5)*eta21 +((1+iss)^6)*eta22 +((1+iss)^7)*eta23 +((1+iss)^8)*eta24   )*(1/iss)*i(+8)  
+  iss*((1+iss)^16)*(eta17 +((1+iss)^1)*eta18 +((1+iss)^2)*eta19 +((1+iss)^3)*eta20 +((1+iss)^4)*eta21  +((1+iss)^5)*eta22 +((1+iss)^6)*eta23 +((1+iss)^7)*eta24    )*(1/iss)*i(+7)  
+  iss*((1+iss)^17)*(eta18 +((1+iss)^1)*eta19 +((1+iss)^2)*eta20 +((1+iss)^3)*eta21 +((1+iss)^4)*eta22  +((1+iss)^5)*eta23 +((1+iss)^6)*eta24    )*(1/iss)*i(+6)  
+  iss*((1+iss)^18)*(eta19 +((1+iss)^1)*eta20 +((1+iss)^2)*eta21 +((1+iss)^3)*eta22 +((1+iss)^4)*eta23  +((1+iss)^5)*eta24    )*(1/iss)*i(+5)  
+  iss*((1+iss)^19)*(eta20 +((1+iss)^1)*eta21 +((1+iss)^2)*eta22 +((1+iss)^3)*eta23 +((1+iss)^4)*eta24    )*(1/iss)*i(+4)  
+  iss*((1+iss)^20)*(eta21 +((1+iss)^1)*eta22 +((1+iss)^2)*eta23 +((1+iss)^3)*eta24    )*(1/iss)*i(+3)  
+  iss*((1+iss)^21)*(eta22 +((1+iss)^1)*eta23 +((1+iss)^2)*eta24    )*(1/iss)*i(+2)  
+  iss*((1+iss)^22)*(eta23 +((1+iss)^1)*eta24    )*(1/iss)*i(+1)  
+  iss*((1+iss)^23)*(eta24     )*(1/iss)*i 
+  eta1*(1+iss)*(p(+23)-p(+24)) + eta2*((1+iss)^2)*(p(+22)-p(+24)) + eta3*((1+iss)^3)*(p(+21)-p(+24)) + eta4*((1+iss)^4)*(p(+20)-p(+24)) + eta5*((1+iss)^5)*(p(+19) -p(+24)) + eta6*((1+iss)^6)*(p(+18) -p(+24)) + eta7*((1+iss)^7)*(p(+17) -p(+24)) + eta8*((1+iss)^8)*(p(+16) -p(+24)) + eta9*((1+iss)^9)*(p(+15) -p(+24)) + eta10*((1+iss)^10)*(p(+14) -p(+24)) + eta11*((1+iss)^11)*(p(+13) -p(+24)) + eta12*((1+iss)^12)*(p(+12) -p(+24)) + eta13*((1+iss)^13)*(p(+11)-p(+24))+ eta14*((1+iss)^14)*(p(+10) -p(+24)) + eta15*((1+iss)^15)*(p(+9) -p(+24))  + eta16*((1+iss)^16)*(p(+8) -p(+24)) + eta17*((1+iss)^17)*(p(+7) -p(+24)) + eta18*((1+iss)^18)*(p(+6) -p(+24)) + eta19*((1+iss)^19)*(p(+5) -p(+24)) + eta20*((1+iss)^20)*(p(+4) -p(+24)) + eta21*((1+iss)^21)*(p(+3) -p(+24)) + eta22*((1+iss)^22)*(p(+2) -p(+24)) + eta23*((1+iss)^23)*(p(+1) -p(+24)) + eta24*((1+iss)^24)*(p -p(+24))  );    

a = rdscale*(1-psi)*nf(-23) + psi*a(-1);
  @#endif
@#endif









@#if ttb_spec==10
 @#if TI==10

 rd=(1+alfa)*(eta1*nf(-9)+eta2*nf(-8)+eta3*nf(-7)+eta4*nf(-6)+eta5*nf(-5)+eta6*nf(-4)+eta7*nf(-3)+eta8*nf(-2)+eta9*nf(-1)+eta10*nf);
 pf(+10)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10)*(1/iss)*i(+9) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10)*(1/iss)*i(+8) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10)*(1/iss)*i(+7) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10)*(1/iss)*i(+6) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10)*(1/iss)*i(+5) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10)*(1/iss)*i(+4) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10)*(1/iss)*i(+3) + iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10)*(1/iss)*i(+2) + iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10)*(1/iss)*i(+1) + iss*((1+iss)^9)*(eta10)*(1/iss)*i + eta1*(1+iss)*(p(+9)-p(+10)) + eta2*((1+iss)^2)*(p(+8)-p(+10)) + eta3*((1+iss)^3)*(p(+7)-p(+10)) + eta4*((1+iss)^4)*(p(+6)-p(+10)) + eta5*((1+iss)^5)*(p(+5) -p(+10)) + eta6*((1+iss)^6)*(p(+4) -p(+10)) + eta7*((1+iss)^7)*(p(+3) -p(+10)) + eta8*((1+iss)^8)*(p(+2) -p(+10)) + eta9*((1+iss)^9)*(p(+1) -p(+10)) + eta10*((1+iss)^10)*(p -p(+10))  );    
 a = rdscale*(1-psi)*nf(-9) + psi*a(-1);
 @#endif
@#endif
 
c-hb*c(-1)= (1+gamma*hb)*(c(+1)-hb*c) - gamma*hb*(c(+2)-hb*c(+1)) - ((1-hb)/sigc)*(1-gamma*hb)*((iss/(1+iss))*(1/iss)*i + p - p(+1) + uc(+1) + ud(+1)) + ((1-hb)/sigc)*(uc + ud + gamma*hb*(uc(+2) +ud(+2))) ; 
sigh*h = w - uh + (1/(1-hb*gamma))*( (-sigc/(1-hb))*(c - hb*c(-1)) + uc ) + (hb*gamma/(1-hb*gamma))*( (sigc/(1-hb))*(c(+1)-hb*c) - uc(+1) - ud(+1) +ud );
sigm*(m-p) = (-1/(1+iss))*(1/iss)*i + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(c-hb*c(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(c(+1)-hb*c) + uc(+1) +ud(+1) - ud) );
lamda = (1-delta)*gamma*lamda(+1) - ((1-(1-delta)*gamma)/(1-hb*gamma))*( (sigc/(1-hb))*(c(+1)-hb*c) - (1/rss)*r(+1) - uc(+1) - ud(+1) - hb*gamma*( (sigc/(1-hb))*(c(+2)-hb*c(+1)) - (1/rss)*r(+1) -uc(+2) -ud(+2) ) );
iv =  (gamma/(1+gamma))*iv(+1) + (1/(1+gamma))*iv(-1) + (1/(ss*(1+gamma)))*lamda  + (1/(ss*(1+gamma)*(1-hb*gamma)))*( (sigc/(1-hb))*(c-hb*c(-1)) - uc - ud +gamma*hb*( (-sigc/(1-hb))*(c(+1)-hb*c) + uc(+1) +ud(+1) ) ) ;
//yss*y = css*c + ivss*iv + rdss*rd;
yss*y = css*c + ivss*iv + rdss*rd +ggss*gg;

(1-delta)*k(-1) +delta*iv - k = 0;
tiv= (ivss*iv + rdss*rd)/(ivss+rdss); //*total investment: R&D investment+ capital investment
y = (1/(phi-1))*a(-1) + u + theta*k(-1) + (1-theta)*h;
h - k(-1) - (1/rss)*r + w = 0;
pai = gamma*(pai(+1)) + (((1-gamma*rho*psi)*(1-rho*psi))/(rho*psi))*(-u + theta*(1/rss)*r +(1-theta)*w - (1/(phi-1))*a(-1)) + (gamma/(phi-1))*(a - a(-1)) - (1/(phi-1))*(a(-1) - a(-2));
pfss*pf - (gamma*psi)*pfss*pf(+1) + (gamma*psi)*pfss*(p - p(+1) + iss*(1/iss)*i) = (1/phi)*yss*y +((phi-1)/phi)*yss*u - ((phi-1)/phi)*(1-theta)*yss*w -  ((phi-1)/phi)*theta*yss*(1/rss)*r;
pai = p - p(-1);
gg=(1/tau)*uf +y; // * govt expenditure

y_g=y-y(-1);
c_g=c-c(-1);
tiv_g=tiv-tiv(-1);
iv_g=iv-iv(-1);
h_g=h-h(-1);
w_g=w-w(-1);
tfp_g=tfp-tfp(-1);

tfpe= 0.999*tfpe(-1)+tfp_g + etfp;

rr=i-pai(+1); //*real interest rate
pf_exp=pf(+8);
///////////////////////////////////////////////////////////////
// Monetary Policy (Inerest rate rule) & interest rate shock //
//////////////////////////////////////////////////////////////
iss*(1/iss)*i = ((1+iss)^(xxx))*(xxx*(iss/(1+iss))*(1/iss)*i(-1) + (1-xxx)*((xpai)*(pai)+(xy)*(yg)) + v);

@#if ishock==1
 v = rhoi*v(-1) + ei;
@#endif

@#if ishock==0
 v = rhoi*v(-1);
@#endif

@#if fshock==1
uf = rhouf*uf(-1)+ef;
@#endif

@#if fshock==0
uf = rhouf*uf(-1)+0.000001*ef;
@#endif

@#if ushock==1
 @#if news==1
   

 @#if eshock==16
    u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9) +eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue(-16);
  
 

 @#if var_data==22 
    conf12= eue(-12); conf11= eue11(-11);conf10= eue10(-10);conf9= eue9(-9);conf8= eue8(-8);
    @#if conf_var==2
    conf=rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4;
   @#endif
 @#endif 

@#endif  
  
   
 @#endif

 @#if news==0

 
 @#if eshock==16
    u = rhou*u(-1); 
  
    @#if var_data==22 
    conf12= 0; conf11=0;conf10= 0;conf9=0;conf8= 0;
   @#if conf_var==2
    conf=0;
   @#endif

  
 @#endif 


   @#endif  
  

 @#endif
@#endif 

@#if ushock==0
 @#if news==1
   
    
   @#if eshock==16
    u = rhou*u(-1) + eue1(-1) + eue2(-2) + eue3(-3) +eue4(-4)+eue5(-5)+eue6(-6)+eue7(-7)+eue8(-8)+eue9(-9) +eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue(-16); 
     

    @#if var_data==22 
    conf12= eue(-12); conf11= eue11(-11);conf10= eue10(-10);conf9= eue9(-9);conf8= eue8(-8);
    
   @#if conf_var==2
     conf=rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4;
   @#endif

 @#endif 


   @#endif  
  @#endif


 @#if news==0



  @#if eshock==16
    u = rhou*u(-1); 
   @#if var_data==22 
    conf12= 0; conf11=0;conf10= 0;conf9=0;conf8= 0;
    @#if conf_var==1
    conf=0;
   @#endif
   @#if conf_var==2
    conf=0;
   @#endif

   @#endif 
  @#endif  
  

 @#endif 
@#endif



 
////////////////////////
// Preference shocks //
//////////////////////
@#if cshock==0
uc = rhouc*uc(-1);
@#endif
@#if cshock==1
uc = rhouc*uc(-1) +euc; 
@#endif

@#if hshock==0
uh = rhouh*uh(-1); //* Preference shock
@#endif
@#if hshock==1
uh = rhouh*uh(-1) +euh; //* Preference shock
@#endif


ud = rhoud*ud(-1);




um = rhoum*um(-1);

yg=y-ynobs; //* gdp gap with measurement errors
hg = h-hn;

pinf= (gamma/(phi-1))*(a - a(-1)) - (1/(phi-1))*(a(-1) - a(-2));
mcgap=(-u + theta*(1/rss)*r +(1-theta)*w - (1/(phi-1))*an(-1)); //* (mc - mcn)
mcgape=mcgap+mee; 
mc=-u + theta*(1/rss)*r +(1-theta)*w; //*real marginal cost under sticky price
//mc=-(1/(phi-1))*a(-1)-u + theta*(1/rss)*r +(1-theta)*w; //*this real marginal cost  under sticky price
mce=mc+mee;
fpai = pai(+1);
ppai = pai(-1);
fpaie = pai(+1)+meee;
paiee=pai+meee+me;//*pai with 2sls errors & measument errors
paie=pai+me;//*pai with measument errors
ppaie=paie(-1);//*ppai with measument errors


inftr =  (((1-gamma*rho*psi)*(1-rho*psi))/(rho*psi))*(- (1/(phi-1))*a(-1)) + (gamma/(phi-1))*(a - a(-1)) - (1/(phi-1))*(a(-1) - a(-2)); //*the residural part of inflation eqn
tfp=(1)*((1/(phi-1))*a(-1) + u) ;
//tfp= u ;

tfpu=u;
shu=eu;
shi=ei;
shc=euc;
shh=euh;
shue1=eue1;
shue3=eue3;
shue4=eue4;
shue8=eue8;
aa=a(-1);

tq=lamda-(1/(1-gamma*hb))*( (sigc/(1-hb))*(-c+hb*c(-1))+uc+ud - gamma*hb*( (sigc/(1-hb))*(-c(+1)+hb*c) + uc(+1)+ud(+1) ) );
sp=tq+k(-1); 
stk=pf+sp;
//*stk is the sum of market value of capital and the intermediate firm *//  
//* q = stock price/capita *//  


//////////////////////////
//* natural variables *//
////////////////////////
 
@#if ttb_spec==8

 
 @#if TI==1
 rdn=(1+alfa)*(eta1*nfn(-7)+eta2*nfn(-6)+eta3*nfn(-5)+eta4*nfn(-4)+eta5*nfn(-3)+eta6*nfn(-2)+eta7*nfn(-1)+eta8*nfn);
 pfn(+8)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8)*(1/iss)*in(+7) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8)*(1/iss)*in(+6) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8)*(1/iss)*in(+5) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8)*(1/iss)*in(+4) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8)*(1/iss)*in(+3) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8)*(1/iss)*in(+2) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8)*(1/iss)*in(+1) + iss*((1+iss)^7)*(eta8)*(1/iss)*in + eta1*(1+iss)*(pn(+7)-pn(+8)) + eta2*((1+iss)^2)*(pn(+6)-pn(+8)) + eta3*((1+iss)^3)*(pn(+5)-pn(+8)) + eta4*((1+iss)^4)*(pn(+4)-pn(+8)) + eta5*((1+iss)^5)*(pn(+3) -pn(+8)) + eta6*((1+iss)^6)*(pn(+2) -pn(+8)) + eta7*((1+iss)^7)*(pn(+1) -pn(+8)) + eta8*((1+iss)^8)*(pn -pn(+8))  );    
 an = rdscale*(1-psi)*nfn(-7) + psi*an(-1);
 @#endif

 @#if TI==0
 rdn=(1+alfa)*(nfn);
 pfn(+1)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* ( (in)  +(1+iss)*(pn-pn(+1)) );    
 an = rdscale*(1-psi)*nfn + psi*an(-1);
 @#endif


@#endif

@#if ttb_spec==12
 @#if TI==12
 rdn=(1+alfa)*(eta1*nfn(-11)+eta2*nfn(-10)+eta3*nfn(-9)+eta4*nfn(-8)+eta5*nfn(-7)+eta6*nfn(-6)+eta7*nfn(-5)+eta8*nfn(-4)+eta9*nfn(-3)+eta10*nfn(-2)+eta11*nfn(-1)+eta12*nfn);
 pfn(+12)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10 +((1+iss)^10)*eta11 +((1+iss)^11)*eta12)*(1/iss)*in(+11) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10 +((1+iss)^9)*eta11 +((1+iss)^10)*eta12)*(1/iss)*in(+10) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10 +((1+iss)^8)*eta11 +((1+iss)^9)*eta12)*(1/iss)*in(+9) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10 +((1+iss)^7)*eta11 +((1+iss)^8)*eta12)*(1/iss)*in(+8) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10 +((1+iss)^6)*eta11 +((1+iss)^7)*eta12)*(1/iss)*in(+7) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10 +((1+iss)^5)*eta11 +((1+iss)^6)*eta12)*(1/iss)*in(+6) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10 +((1+iss)^4)*eta11+((1+iss)^5)*eta12)*(1/iss)*in(+5) + iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10 +((1+iss)^3)*eta11 +((1+iss)^4)*eta12)*(1/iss)*in(+4) + iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10 +((1+iss)^2)*eta11 +((1+iss)^3)*eta12)*(1/iss)*in(+3) + iss*((1+iss)^9)*(eta10 +((1+iss)^1)*eta11 +((1+iss)^2)*eta12)*(1/iss)*in(+2)  + iss*((1+iss)^10)*(eta11 +((1+iss)^1)*eta12)*(1/iss)*in(+1)  + iss*((1+iss)^11)*(eta12)*(1/iss)*in + eta1*(1+iss)*(pn(+11)-pn(+12)) + eta2*((1+iss)^2)*(pn(+10)-pn(+12)) + eta3*((1+iss)^3)*(pn(+9)-pn(+12)) + eta4*((1+iss)^4)*(pn(+8)-pn(+12)) + eta5*((1+iss)^5)*(pn(+7) -pn(+12)) + eta6*((1+iss)^6)*(pn(+6) -pn(+12)) + eta7*((1+iss)^7)*(pn(+5) -pn(+12)) + eta8*((1+iss)^8)*(pn(+4) -pn(+12)) + eta9*((1+iss)^9)*(pn(+3) -pn(+12)) + eta10*((1+iss)^10)*(pn(+2) -pn(+12)) + eta11*((1+iss)^11)*(pn(+1) -pn(+12)) + eta12*((1+iss)^12)*(pn -pn(+12))  );    
 an = rdscale*(1-psi)*nfn(-11) + psi*an(-1);
 @#endif
@#endif


@#if ttb_spec==16
 @#if TI==16
rdn=(1+alfa)*(eta1*nfn(-15)+eta2*nfn(-14)+eta3*nfn(-13)+eta4*nfn(-12)+eta5*nfn(-11)+eta6*nfn(-10)+eta7*nfn(-9)+eta8*nfn(-8)+eta9*nfn(-7)+eta10*nfn(-6)+eta11*nfn(-5)+eta12*nfn*(-4)+eta13*nfn*(-3)+eta14*nfn*(-2)+eta15*nfn*(-1)+eta16*nfn);
 pfn(+16)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* 
(  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10 +((1+iss)^10)*eta11 +((1+iss)^11)*eta12+((1+iss)^12)*eta13+((1+iss)^13)*eta14+((1+iss)^14)*eta15+((1+iss)^15)*eta16)*(1/iss)*in(+15) 
+  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10 +((1+iss)^9)*eta11 +((1+iss)^10)*eta12 +((1+iss)^11)*eta13 +((1+iss)^12)*eta14 +((1+iss)^13)*eta15 +((1+iss)^14)*eta16)*(1/iss)*in(+14) 
+  iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10 +((1+iss)^8)*eta11 +((1+iss)^9)*eta12 +((1+iss)^10)*eta13  +((1+iss)^11)*eta14 +((1+iss)^12)*eta15 +((1+iss)^13)*eta16)*(1/iss)*in(+13) 
+  iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10 +((1+iss)^7)*eta11 +((1+iss)^8)*eta12 +((1+iss)^9)*eta13 +((1+iss)^10)*eta14 +((1+iss)^11)*eta15 +((1+iss)^12)*eta16 )*(1/iss)*in(+12) 
+  iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10 +((1+iss)^6)*eta11 +((1+iss)^7)*eta12 +((1+iss)^8)*eta13 +((1+iss)^9)*eta14 +((1+iss)^10)*eta15 +((1+iss)^11)*eta16)*(1/iss)*in(+11) 
+  iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10 +((1+iss)^5)*eta11 +((1+iss)^6)*eta12 +((1+iss)^7)*eta13 +((1+iss)^8)*eta14 +((1+iss)^9)*eta15 +((1+iss)^10)*eta16)*(1/iss)*in(+10) 
+  iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10 +((1+iss)^4)*eta11+((1+iss)^5)*eta12 +((1+iss)^6)*eta13 +((1+iss)^7)*eta14 +((1+iss)^8)*eta15 +((1+iss)^9)*eta16)*(1/iss)*in(+9) 
+  iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10 +((1+iss)^3)*eta11 +((1+iss)^4)*eta12  +((1+iss)^5)*eta13  +((1+iss)^6)*eta14  +((1+iss)^7)*eta15  +((1+iss)^8)*eta16)*(1/iss)*in(+8) 
+  iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10 +((1+iss)^2)*eta11 +((1+iss)^3)*eta12 +((1+iss)^4)*eta13  +((1+iss)^5)*eta14  +((1+iss)^6)*eta15  +((1+iss)^7)*eta16 )*(1/iss)*in(+7) 
+  iss*((1+iss)^9)*(eta10 +((1+iss)^1)*eta11 +((1+iss)^2)*eta12 +((1+iss)^3)*eta13 +((1+iss)^4)*eta14  +((1+iss)^5)*eta15  +((1+iss)^6)*eta16  )*(1/iss)*in(+6) 
+  iss*((1+iss)^10)*(eta11 +((1+iss)^1)*eta12 +((1+iss)^2)*eta13 +((1+iss)^3)*eta14 +((1+iss)^4)*eta15  +((1+iss)^5)*eta16    )*(1/iss)*in(+5) 
+  iss*((1+iss)^11)*(eta12 +((1+iss)^1)*eta13 +((1+iss)^2)*eta14 +((1+iss)^3)*eta15 +((1+iss)^4)*eta16      )*(1/iss)*in(+4) 
+  iss*((1+iss)^12)*(eta13 +((1+iss)^1)*eta14 +((1+iss)^2)*eta15 +((1+iss)^3)*eta16    )*(1/iss)*in(+3) 
+  iss*((1+iss)^13)*(eta14 +((1+iss)^1)*eta15 +((1+iss)^2)*eta16    )*(1/iss)*in(+2)  
+  iss*((1+iss)^14)*(eta15 +((1+iss)^1)*eta16    )*(1/iss)*in(+1)  
+  iss*((1+iss)^15)*(eta16   )*(1/iss)*in 
+  eta1*(1+iss)*(pn(+15)-pn(+16)) + eta2*((1+iss)^2)*(pn(+14)-pn(+16)) + eta3*((1+iss)^3)*(pn(+13)-pn(+16)) + eta4*((1+iss)^4)*(pn(+12)-pn(+16)) + eta5*((1+iss)^5)*(pn(+11) -pn(+16)) + eta6*((1+iss)^6)*(pn(+10) -pn(+16)) + eta7*((1+iss)^7)*(pn(+9) -pn(+16)) + eta8*((1+iss)^8)*(pn(+8) -pn(+16)) + eta9*((1+iss)^9)*(pn(+7) -pn(+16)) + eta10*((1+iss)^10)*(pn(+6) -pn(+16)) + eta11*((1+iss)^11)*(pn(+5) -pn(+16)) + eta12*((1+iss)^12)*(pn(+4) -pn(+16)) + eta13*((1+iss)^13)*(pn(+3) -pn(+16))+ eta14*((1+iss)^14)*(pn(+2) -pn(+16))+ eta15*((1+iss)^15)*(pn(+1) -pn(+16))+ eta16*((1+iss)^16)*(pn -pn(+16))  );    
an = rdscale*(1-psi)*nfn(-15) + psi*an(-1);
@#endif
@#endif


@#if ttb_spec==20
 @#if TI==20

 rdn=(1+alfa)*(eta1*nfn(-19)+eta2*nfn(-18)+eta3*nfn(-17)+eta4*nfn(-16)+eta5*nfn(-15)+eta6*nfn(-14)+eta7*nfn(-13)+eta8*nfn(-12)+eta9*nfn(-11)+eta10*nfn(-10)+eta11*nfn(-9)+eta12*nfn(-8)+eta13*nfn(-7)+eta14*nfn(-6)+eta15*nfn(-5)+eta16*nfn*(-4)+eta17*nfn*(-3)+eta18*nfn*(-2)+eta19*nfn*(-1)+eta20*nfn);          
 pfn(+20)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* 
(  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10 +((1+iss)^10)*eta11 +((1+iss)^11)*eta12+((1+iss)^12)*eta13+((1+iss)^13)*eta14 +((1+iss)^14)*eta15+((1+iss)^15)*eta16+((1+iss)^16)*eta17+((1+iss)^17)*eta18+((1+iss)^18)*eta19 +((1+iss)^19)*eta20  )*(1/iss)*in(+19) 
+  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10 +((1+iss)^9)*eta11 +((1+iss)^10)*eta12 +((1+iss)^11)*eta13 +((1+iss)^12)*eta14 +((1+iss)^13)*eta15 +((1+iss)^14)*eta16  +((1+iss)^15)*eta17 +((1+iss)^16)*eta18 +((1+iss)^17)*eta19 +((1+iss)^18)*eta20  )*(1/iss)*in(+18) 
+  iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10 +((1+iss)^8)*eta11 +((1+iss)^9)*eta12 +((1+iss)^10)*eta13  +((1+iss)^11)*eta14 +((1+iss)^12)*eta15 +((1+iss)^13)*eta16 +((1+iss)^14)*eta17 +((1+iss)^15)*eta18 +((1+iss)^16)*eta19 +((1+iss)^17)*eta20  )*(1/iss)*in(+17) 
+  iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10 +((1+iss)^7)*eta11 +((1+iss)^8)*eta12 +((1+iss)^9)*eta13 +((1+iss)^10)*eta14 +((1+iss)^11)*eta15 +((1+iss)^12)*eta16 +((1+iss)^13)*eta17 +((1+iss)^14)*eta18 +((1+iss)^15)*eta19 +((1+iss)^16)*eta20  )*(1/iss)*in(+16) 
+  iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10 +((1+iss)^6)*eta11 +((1+iss)^7)*eta12 +((1+iss)^8)*eta13 +((1+iss)^9)*eta14 +((1+iss)^10)*eta15 +((1+iss)^11)*eta16 +((1+iss)^12)*eta17 +((1+iss)^13)*eta18 +((1+iss)^14)*eta19 +((1+iss)^15)*eta20  )*(1/iss)*in(+15) 

+  iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10 +((1+iss)^5)*eta11 +((1+iss)^6)*eta12 +((1+iss)^7)*eta13 +((1+iss)^8)*eta14 +((1+iss)^9)*eta15 +((1+iss)^10)*eta16 +((1+iss)^11)*eta17 +((1+iss)^12)*eta18 +((1+iss)^13)*eta19 +((1+iss)^14)*eta20  )*(1/iss)*in(+14) 
+  iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10 +((1+iss)^4)*eta11+((1+iss)^5)*eta12 +((1+iss)^6)*eta13 +((1+iss)^7)*eta14 +((1+iss)^8)*eta15   +((1+iss)^9)*eta16 +((1+iss)^10)*eta17 +((1+iss)^11)*eta18 +((1+iss)^12)*eta19 +((1+iss)^13)*eta20  )*(1/iss)*in(+13) 
+  iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10 +((1+iss)^3)*eta11 +((1+iss)^4)*eta12  +((1+iss)^5)*eta13  +((1+iss)^6)*eta14  +((1+iss)^7)*eta15 +((1+iss)^8)*eta16 +((1+iss)^9)*eta17 +((1+iss)^10)*eta18 +((1+iss)^11)*eta19 +((1+iss)^12)*eta20  )*(1/iss)*in(+12) 
+  iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10 +((1+iss)^2)*eta11 +((1+iss)^3)*eta12 +((1+iss)^4)*eta13  +((1+iss)^5)*eta14  +((1+iss)^6)*eta15  +((1+iss)^7)*eta16 +((1+iss)^8)*eta17 +((1+iss)^9)*eta18 +((1+iss)^10)*eta19 +((1+iss)^11)*eta20   )*(1/iss)*in(+11) 
+  iss*((1+iss)^9)*(eta10 +((1+iss)^1)*eta11 +((1+iss)^2)*eta12 +((1+iss)^3)*eta13 +((1+iss)^4)*eta14  +((1+iss)^5)*eta15  +((1+iss)^6)*eta16 +((1+iss)^7)*eta17 +((1+iss)^8)*eta18 +((1+iss)^9)*eta19 +((1+iss)^10)*eta20    )*(1/iss)*in(+10) 
+  iss*((1+iss)^10)*(eta11 +((1+iss)^1)*eta12 +((1+iss)^2)*eta13 +((1+iss)^3)*eta14 +((1+iss)^4)*eta15  +((1+iss)^5)*eta16 +((1+iss)^6)*eta17 +((1+iss)^7)*eta18 +((1+iss)^8)*eta19 +((1+iss)^9)*eta20     )*(1/iss)*in(+9) 
+  iss*((1+iss)^11)*(eta12 +((1+iss)^1)*eta13 +((1+iss)^2)*eta14 +((1+iss)^3)*eta15 +((1+iss)^4)*eta16 +((1+iss)^5)*eta17 +((1+iss)^6)*eta18 +((1+iss)^7)*eta19 +((1+iss)^8)*eta20     )*(1/iss)*in(+8) 
+  iss*((1+iss)^12)*(eta13 +((1+iss)^1)*eta14 +((1+iss)^2)*eta15 +((1+iss)^3)*eta16 +((1+iss)^4)*eta17 +((1+iss)^5)*eta18 +((1+iss)^6)*eta19 +((1+iss)^7)*eta20    )*(1/iss)*in(+7) 
+  iss*((1+iss)^13)*(eta14 +((1+iss)^1)*eta15 +((1+iss)^2)*eta16 +((1+iss)^3)*eta17 +((1+iss)^4)*eta18 +((1+iss)^5)*eta19 +((1+iss)^6)*eta20    )*(1/iss)*in(+6)  
+  iss*((1+iss)^14)*(eta15 +((1+iss)^1)*eta16 +((1+iss)^2)*eta17 +((1+iss)^3)*eta18 +((1+iss)^4)*eta19 +((1+iss)^5)*eta20   )*(1/iss)*in(+5)  
+  iss*((1+iss)^15)*(eta16 +((1+iss)^1)*eta17 +((1+iss)^2)*eta18 +((1+iss)^3)*eta19 +((1+iss)^4)*eta20  )*(1/iss)*in(+4)  
+  iss*((1+iss)^16)*(eta17 +((1+iss)^1)*eta18 +((1+iss)^2)*eta19 +((1+iss)^3)*eta20    )*(1/iss)*in(+3)  
+  iss*((1+iss)^17)*(eta18 +((1+iss)^1)*eta19 +((1+iss)^2)*eta20   )*(1/iss)*in(+2)  
+  iss*((1+iss)^18)*(eta19 +((1+iss)^1)*eta20   )*(1/iss)*in(+1)  
+  iss*((1+iss)^19)*(eta20     )*(1/iss)*in 
+  eta1*(1+iss)*(pn(+19)-pn(+20)) + eta2*((1+iss)^2)*(pn(+18)-pn(+20)) + eta3*((1+iss)^3)*(pn(+17)-pn(+20)) + eta4*((1+iss)^4)*(pn(+16)-pn(+20)) + eta5*((1+iss)^5)*(pn(+15) -pn(+20)) + eta6*((1+iss)^6)*(pn(+14) -pn(+20)) + eta7*((1+iss)^7)*(pn(+13) -pn(+20)) + eta8*((1+iss)^8)*(pn(+12) -pn(+20)) + eta9*((1+iss)^9)*(pn(+11) -pn(+20)) + eta10*((1+iss)^10)*(pn(+10) -pn(+20)) + eta11*((1+iss)^11)*(pn(+9) -pn(+20)) + eta12*((1+iss)^12)*(pn(+8) -pn(+20)) + eta13*((1+iss)^13)*(pn(+7)-pn(+20))+ eta14*((1+iss)^14)*(pn(+6) -pn(+20)) + eta15*((1+iss)^15)*(pn(+5) -pn(+20))  + eta16*((1+iss)^16)*(pn(+4) -pn(+20)) + eta17*((1+iss)^17)*(pn(+3) -pn(+20)) + eta18*((1+iss)^18)*(pn(+2) -pn(+20)) + eta19*((1+iss)^19)*(pn(+1) -pn(+20)) + eta20*((1+iss)^20)*(pn -pn(+20))  );    

an = rdscale*(1-psi)*nfn(-19) + psi*an(-1);
  @#endif
@#endif






@#if ttb_spec==24
 @#if TI==24
rdn=(1+alfa)*(eta1*nfn(-23)+eta2*nfn(-22)+eta3*nfn(-21)+eta4*nfn(-20)+eta5*nfn(-19)+eta6*nfn(-18)+eta7*nfn(-17)+eta8*nfn(-16)+eta9*nfn(-15)+eta10*nfn(-14)+eta11*nfn(-13)+eta12*nfn(-12)+eta13*nfn(-11)+eta14*nf(-10)+eta15*nfn(-9)+eta16*nfn(-8)+eta17*nfn(-7)+eta18*nfn(-6)+eta19*nfn(-5)+eta20*nfn*(-4)+eta21*nfn*(-3)+eta22*nfn*(-2)+eta23*nfn*(-1)+eta24*nfn);          
 pfn(+24)=alfa*nf+((nfss^alfa)/(epsilon*x*pfss))* 
(  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10 +((1+iss)^10)*eta11 +((1+iss)^11)*eta12+((1+iss)^12)*eta13+((1+iss)^13)*eta14 +((1+iss)^14)*eta15+((1+iss)^15)*eta16+((1+iss)^16)*eta17+((1+iss)^17)*eta18+((1+iss)^18)*eta19 +((1+iss)^19)*eta20  +((1+iss)^20)*eta21 +((1+iss)^21)*eta22 +((1+iss)^22)*eta23 +((1+iss)^23)*eta24)*(1/iss)*in(+23) 
+  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10 +((1+iss)^9)*eta11 +((1+iss)^10)*eta12 +((1+iss)^11)*eta13 +((1+iss)^12)*eta14 +((1+iss)^13)*eta15 +((1+iss)^14)*eta16  +((1+iss)^15)*eta17 +((1+iss)^16)*eta18 +((1+iss)^17)*eta19 +((1+iss)^18)*eta20  +((1+iss)^19)*eta21 +((1+iss)^20)*eta22 +((1+iss)^21)*eta23 +((1+iss)^22)*eta24)*(1/iss)*in(+22) 
+  iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10 +((1+iss)^8)*eta11 +((1+iss)^9)*eta12 +((1+iss)^10)*eta13  +((1+iss)^11)*eta14 +((1+iss)^12)*eta15 +((1+iss)^13)*eta16 +((1+iss)^14)*eta17 +((1+iss)^15)*eta18 +((1+iss)^16)*eta19 +((1+iss)^17)*eta20  +((1+iss)^18)*eta21 +((1+iss)^19)*eta22 +((1+iss)^20)*eta23 +((1+iss)^21)*eta24)*(1/iss)*in(+21) 
+  iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10 +((1+iss)^7)*eta11 +((1+iss)^8)*eta12 +((1+iss)^9)*eta13 +((1+iss)^10)*eta14 +((1+iss)^11)*eta15 +((1+iss)^12)*eta16 +((1+iss)^13)*eta17 +((1+iss)^14)*eta18 +((1+iss)^15)*eta19 +((1+iss)^16)*eta20  +((1+iss)^17)*eta21 +((1+iss)^18)*eta22 +((1+iss)^19)*eta23 +((1+iss)^20)*eta24)*(1/iss)*in(+20) 
+  iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10 +((1+iss)^6)*eta11 +((1+iss)^7)*eta12 +((1+iss)^8)*eta13 +((1+iss)^9)*eta14 +((1+iss)^10)*eta15 +((1+iss)^11)*eta16 +((1+iss)^12)*eta17 +((1+iss)^13)*eta18 +((1+iss)^14)*eta19 +((1+iss)^15)*eta20  +((1+iss)^16)*eta21 +((1+iss)^17)*eta22 +((1+iss)^18)*eta23 +((1+iss)^19)*eta24)*(1/iss)*in(+19) 

+  iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10 +((1+iss)^5)*eta11 +((1+iss)^6)*eta12 +((1+iss)^7)*eta13 +((1+iss)^8)*eta14 +((1+iss)^9)*eta15 +((1+iss)^10)*eta16 +((1+iss)^11)*eta17 +((1+iss)^12)*eta18 +((1+iss)^13)*eta19 +((1+iss)^14)*eta20  +((1+iss)^15)*eta21 +((1+iss)^16)*eta22 +((1+iss)^17)*eta23 +((1+iss)^18)*eta24)*(1/iss)*in(+18) 
+  iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10 +((1+iss)^4)*eta11+((1+iss)^5)*eta12 +((1+iss)^6)*eta13 +((1+iss)^7)*eta14 +((1+iss)^8)*eta15   +((1+iss)^9)*eta16 +((1+iss)^10)*eta17 +((1+iss)^11)*eta18 +((1+iss)^12)*eta19 +((1+iss)^13)*eta20  +((1+iss)^14)*eta21 +((1+iss)^15)*eta22 +((1+iss)^16)*eta23 +((1+iss)^17)*eta24)*(1/iss)*in(+17) 
+  iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10 +((1+iss)^3)*eta11 +((1+iss)^4)*eta12  +((1+iss)^5)*eta13  +((1+iss)^6)*eta14  +((1+iss)^7)*eta15 +((1+iss)^8)*eta16 +((1+iss)^9)*eta17 +((1+iss)^10)*eta18 +((1+iss)^11)*eta19 +((1+iss)^12)*eta20  +((1+iss)^13)*eta21 +((1+iss)^14)*eta22 +((1+iss)^15)*eta23 +((1+iss)^16)*eta24)*(1/iss)*in(+16) 
+  iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10 +((1+iss)^2)*eta11 +((1+iss)^3)*eta12 +((1+iss)^4)*eta13  +((1+iss)^5)*eta14  +((1+iss)^6)*eta15  +((1+iss)^7)*eta16 +((1+iss)^8)*eta17 +((1+iss)^9)*eta18 +((1+iss)^10)*eta19 +((1+iss)^11)*eta20  +((1+iss)^12)*eta21 +((1+iss)^13)*eta22 +((1+iss)^14)*eta23 +((1+iss)^15)*eta24 )*(1/iss)*in(+15) 
+  iss*((1+iss)^9)*(eta10 +((1+iss)^1)*eta11 +((1+iss)^2)*eta12 +((1+iss)^3)*eta13 +((1+iss)^4)*eta14  +((1+iss)^5)*eta15  +((1+iss)^6)*eta16 +((1+iss)^7)*eta17 +((1+iss)^8)*eta18 +((1+iss)^9)*eta19 +((1+iss)^10)*eta20  +((1+iss)^11)*eta21 +((1+iss)^12)*eta22 +((1+iss)^13)*eta23 +((1+iss)^14)*eta24  )*(1/iss)*in(+14) 
+  iss*((1+iss)^10)*(eta11 +((1+iss)^1)*eta12 +((1+iss)^2)*eta13 +((1+iss)^3)*eta14 +((1+iss)^4)*eta15  +((1+iss)^5)*eta16 +((1+iss)^6)*eta17 +((1+iss)^7)*eta18 +((1+iss)^8)*eta19 +((1+iss)^9)*eta20  +((1+iss)^10)*eta21 +((1+iss)^11)*eta22 +((1+iss)^12)*eta23 +((1+iss)^13)*eta24    )*(1/iss)*in(+13) 
+  iss*((1+iss)^11)*(eta12 +((1+iss)^1)*eta13 +((1+iss)^2)*eta14 +((1+iss)^3)*eta15 +((1+iss)^4)*eta16 +((1+iss)^5)*eta17 +((1+iss)^6)*eta18 +((1+iss)^7)*eta19 +((1+iss)^8)*eta20  +((1+iss)^9)*eta21 +((1+iss)^10)*eta22 +((1+iss)^11)*eta23 +((1+iss)^12)*eta24      )*(1/iss)*in(+12) 
+  iss*((1+iss)^12)*(eta13 +((1+iss)^1)*eta14 +((1+iss)^2)*eta15 +((1+iss)^3)*eta16 +((1+iss)^4)*eta17 +((1+iss)^5)*eta18 +((1+iss)^6)*eta19 +((1+iss)^7)*eta20  +((1+iss)^8)*eta21 +((1+iss)^9)*eta22 +((1+iss)^10)*eta23 +((1+iss)^11)*eta24    )*(1/iss)*in(+11) 
+  iss*((1+iss)^13)*(eta14 +((1+iss)^1)*eta15 +((1+iss)^2)*eta16 +((1+iss)^3)*eta17 +((1+iss)^4)*eta18 +((1+iss)^5)*eta19 +((1+iss)^6)*eta20  +((1+iss)^7)*eta21 +((1+iss)^8)*eta22 +((1+iss)^9)*eta23 +((1+iss)^10)*eta24    )*(1/iss)*in(+10)  
+  iss*((1+iss)^14)*(eta15 +((1+iss)^1)*eta16 +((1+iss)^2)*eta17 +((1+iss)^3)*eta18 +((1+iss)^4)*eta19 +((1+iss)^5)*eta20  +((1+iss)^6)*eta21 +((1+iss)^7)*eta22 +((1+iss)^8)*eta23 +((1+iss)^9)*eta24    )*(1/iss)*in(+9)  
+  iss*((1+iss)^15)*(eta16 +((1+iss)^1)*eta17 +((1+iss)^2)*eta18 +((1+iss)^3)*eta19 +((1+iss)^4)*eta20  +((1+iss)^5)*eta21 +((1+iss)^6)*eta22 +((1+iss)^7)*eta23 +((1+iss)^8)*eta24   )*(1/iss)*in(+8)  
+  iss*((1+iss)^16)*(eta17 +((1+iss)^1)*eta18 +((1+iss)^2)*eta19 +((1+iss)^3)*eta20 +((1+iss)^4)*eta21  +((1+iss)^5)*eta22 +((1+iss)^6)*eta23 +((1+iss)^7)*eta24    )*(1/iss)*in(+7)  
+  iss*((1+iss)^17)*(eta18 +((1+iss)^1)*eta19 +((1+iss)^2)*eta20 +((1+iss)^3)*eta21 +((1+iss)^4)*eta22  +((1+iss)^5)*eta23 +((1+iss)^6)*eta24    )*(1/iss)*in(+6)  
+  iss*((1+iss)^18)*(eta19 +((1+iss)^1)*eta20 +((1+iss)^2)*eta21 +((1+iss)^3)*eta22 +((1+iss)^4)*eta23  +((1+iss)^5)*eta24    )*(1/iss)*in(+5)  
+  iss*((1+iss)^19)*(eta20 +((1+iss)^1)*eta21 +((1+iss)^2)*eta22 +((1+iss)^3)*eta23 +((1+iss)^4)*eta24    )*(1/iss)*in(+4)  
+  iss*((1+iss)^20)*(eta21 +((1+iss)^1)*eta22 +((1+iss)^2)*eta23 +((1+iss)^3)*eta24    )*(1/iss)*in(+3)  
+  iss*((1+iss)^21)*(eta22 +((1+iss)^1)*eta23 +((1+iss)^2)*eta24    )*(1/iss)*in(+2)  
+  iss*((1+iss)^22)*(eta23 +((1+iss)^1)*eta24    )*(1/iss)*in(+1)  
+  iss*((1+iss)^23)*(eta24     )*(1/iss)*in 
+  eta1*(1+iss)*(pn(+23)-pn(+24)) + eta2*((1+iss)^2)*(pn(+22)-pn(+24)) + eta3*((1+iss)^3)*(pn(+21)-pn(+24)) + eta4*((1+iss)^4)*(pn(+20)-pn(+24)) + eta5*((1+iss)^5)*(pn(+19) -pn(+24)) + eta6*((1+iss)^6)*(pn(+18) -pn(+24)) + eta7*((1+iss)^7)*(pn(+17) -pn(+24)) + eta8*((1+iss)^8)*(pn(+16) -pn(+24)) + eta9*((1+iss)^9)*(pn(+15) -pn(+24)) + eta10*((1+iss)^10)*(pn(+14) -pn(+24)) + eta11*((1+iss)^11)*(pn(+13) -pn(+24)) + eta12*((1+iss)^12)*(pn(+12) -pn(+24)) + eta13*((1+iss)^13)*(pn(+11)-pn(+24))+ eta14*((1+iss)^14)*(pn(+10) -pn(+24)) + eta15*((1+iss)^15)*(pn(+9) -pn(+24))  + eta16*((1+iss)^16)*(pn(+8) -pn(+24)) + eta17*((1+iss)^17)*(pn(+7) -pn(+24)) + eta18*((1+iss)^18)*(pn(+6) -pn(+24)) + eta19*((1+iss)^19)*(pn(+5) -pn(+24)) + eta20*((1+iss)^20)*(pn(+4) -pn(+24)) + eta21*((1+iss)^21)*(pn(+3) -pn(+24)) + eta22*((1+iss)^22)*(pn(+2) -pn(+24)) + eta23*((1+iss)^23)*(pn(+1) -pn(+24)) + eta24*((1+iss)^24)*(pn -pn(+24))  );    

an = rdscale*(1-psi)*nfn(-23) + psi*an(-1);
@#endif
@#endif



@#if ttb_spec==10
 @#if TI==10
 rdn=(1+alfa)*(eta1*nfn(-9)+eta2*nfn(-8)+eta3*nfn(-7)+eta4*nfn(-6)+eta5*nfn(-5)+eta6*nfn(-4)+eta7*nfn(-3)+eta8*nfn(-2)+eta9*nfn(-1)+eta10*nfn);
 pfn(+10)=alfa*nfn+((nfss^alfa)/(epsilon*x*pfss))* (  iss*(eta1+(1+iss)*eta2+((1+iss)^2)*eta3+((1+iss)^3)*eta4+((1+iss)^4)*eta5+((1+iss)^5)*eta6 +((1+iss)^6)*eta7 +((1+iss)^7)*eta8 +((1+iss)^8)*eta9 +((1+iss)^9)*eta10)*(1/iss)*in(+9) +  iss*(1+iss)*(eta2+((1+iss)^1)*eta3+((1+iss)^2)*eta4+((1+iss)^3)*eta5+((1+iss)^4)*eta6 +((1+iss)^5)*eta7 +((1+iss)^6)*eta8 +((1+iss)^7)*eta9 +((1+iss)^8)*eta10)*(1/iss)*in(+8) + iss*((1+iss)^2)*(eta3+((1+iss)^1)*eta4+((1+iss)^2)*eta5 +((1+iss)^3)*eta6 +((1+iss)^4)*eta7 +((1+iss)^5)*eta8 +((1+iss)^6)*eta9 +((1+iss)^7)*eta10)*(1/iss)*in(+7) + iss*((1+iss)^3)*(eta4+((1+iss)^1)*eta5 +((1+iss)^2)*eta6 +((1+iss)^3)*eta7 +((1+iss)^4)*eta8 +((1+iss)^5)*eta9 +((1+iss)^6)*eta10)*(1/iss)*in(+6) + iss*((1+iss)^4)*(eta5 +((1+iss)^1)*eta6 +((1+iss)^2)*eta7 +((1+iss)^3)*eta8 +((1+iss)^4)*eta9 +((1+iss)^5)*eta10)*(1/iss)*in(+5) + iss*((1+iss)^5)*(eta6 +((1+iss)^1)*eta7 +((1+iss)^2)*eta8 +((1+iss)^3)*eta9 +((1+iss)^4)*eta10)*(1/iss)*in(+4) + iss*((1+iss)^6)*(eta7 +((1+iss)^1)*eta8 +((1+iss)^2)*eta9 +((1+iss)^3)*eta10)*(1/iss)*in(+3) + iss*((1+iss)^7)*(eta8 +((1+iss)^1)*eta9 +((1+iss)^2)*eta10)*(1/iss)*in(+2) + iss*((1+iss)^8)*(eta9 +((1+iss)^1)*eta10)*(1/iss)*in(+1) + iss*((1+iss)^9)*(eta10)*(1/iss)*in + eta1*(1+iss)*(pn(+9)-pn(+10)) + eta2*((1+iss)^2)*(pn(+8)-pn(+10)) + eta3*((1+iss)^3)*(pn(+7)-pn(+10)) + eta4*((1+iss)^4)*(pn(+6)-pn(+10)) + eta5*((1+iss)^5)*(pn(+5) -pn(+10)) + eta6*((1+iss)^6)*(pn(+4) -pn(+10)) + eta7*((1+iss)^7)*(pn(+3) -pn(+10)) + eta8*((1+iss)^8)*(pn(+2) -pn(+10)) + eta9*((1+iss)^9)*(pn(+1) -pn(+10)) + eta10*((1+iss)^10)*(pn -pn(+10))  );    
 an = rdscale*(1-psi)*nfn(-9) + psi*an(-1);
 @#endif
@#endif
 
cn-hb*cn(-1)= (1+gamma*hb)*(cn(+1)-hb*cn) - gamma*hb*(cn(+2)-hb*cn(+1)) - ((1-hb)/sigc)*(1-gamma*hb)*((iss/(1+iss))*(1/iss)*in + pn - pn(+1) + uc(+1) + ud(+1)) + ((1-hb)/sigc)*(uc + ud + gamma*hb*(uc(+2) +ud(+2))) ; 
sigh*hn = wn - uh + (1/(1-hb*gamma))*( (-sigc/(1-hb))*(cn - hb*cn(-1)) + uc ) + (hb*gamma/(1-hb*gamma))*( (sigc/(1-hb))*(cn(+1)-hb*cn) - uc(+1) - ud(+1) +ud );
sigm*(0-pn) = (-1/(1+iss))*(1/iss)*in + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(cn-hb*cn(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) - ud) );
lamdan = (1-delta)*gamma*lamdan(+1) - ((1-(1-delta)*gamma)/(1-hb*gamma))*( (sigc/(1-hb))*(cn(+1)-hb*cn) - (1/rss)*rn(+1) - uc(+1) - ud(+1) - hb*gamma*( (sigc/(1-hb))*(cn(+2)-hb*cn(+1)) - (1/rss)*rn(+1) -uc(+2) -ud(+2) ) );
ivn =  (gamma/(1+gamma))*ivn(+1) + (1/(1+gamma))*ivn(-1) + (1/(ss*(1+gamma)))*lamdan  + (1/(ss*(1+gamma)*(1-hb*gamma)))*( (sigc/(1-hb))*(cn-hb*cn(-1)) - uc - ud +gamma*hb*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) ) ) ;
//yss*yn = css*cn + ivss*ivn + rdss*rdn;
yss*yn = css*cn + ivss*ivn + rdss*rdn +ggss*ggn;

(1-delta)*kn(-1) +delta*ivn - kn = 0;
tivn= (ivss*ivn + rdss*rdn)/(ivss+rdss);
ynobs=yn;
hn - kn(-1) - (1/rss)*rn + wn = 0;
yn = (1/(phi-1))*an(-1) + u + theta*kn(-1) + (1-theta)*hn;
0 = (((1-gamma*0*psi)*(1-0*psi)))*(-u + theta*(1/rss)*rn +(1-theta)*wn - (1/(phi-1))*an(-1)) ;
pfss*pfn - (gamma*psi)*pfss*pfn(+1) + (gamma*psi)*pfss*(pn - pn(+1) + iss*(1/iss)*in) = (1/phi)*yss*yn +((phi-1)/phi)*yss*u - ((phi-1)/phi)*(1-theta)*yss*wn -  ((phi-1)/phi)*theta*yss*(1/rss)*rn;
pain = pn - pn(-1);

tfpn=((1/(phi-1))*an(-1) + u);
aan=an(-1);
tqn=lamdan-(1/(1-gamma*hb))*( (sigc/(1-hb))*(-cn+hb*cn(-1))+uc+ud - gamma*hb*( (sigc/(1-hb))*(-cn(+1)+hb*cn) + uc(+1)+ud(+1) ) );
//tqn=-(lamdan-(1/(1-gamma*hb))*( (sigc/(1-hb))*(-cn+hb*cn(-1))+uc+ud - gamma*hb*( (sigc/(1-hb))*(-cn(+1)+hb*cn) + uc(+1)+ud(+1) ) ));

mcn=-u + theta*(1/rss)*rn +(1-theta)*wn;
ggn= (1/tau)*uf +yn;

end;

initval;
ynobs=0;
tiv=0;
p = 0; 
pai = 0;
u = 0;
ud= 0;
uc= 0;
uh= 0;
um= 0;
gg=0;
uf=0;

a = 0;
i=0;
lamda =0;
iv=0;

v=0;


nf=0;
nfn=0;


y =0 ;
k = 0;
h = 0;
iv=0;
c = 0;
m = 0;
w =0 ;
r = 0;
rd = 0;
pf =0;
//q=0;
//eue=0;

yn =0 ;
kn = 0;
hn = 0;
cn = 0;
ggn=0;
//mn = 0;
wn =0 ;
rn = 0;
rdn = 0;
pfn =0;
pn = 0; 
pain = 0;
an = 0;
in=0;
ivn=0; 
lamdan=0;
//qn=0;


hg=0;
yg=0;

end;

steady_state_model;
ynobs=0;
tiv=0;
p = 0; 
pai = 0;
u = 0;
ud= 0;
uc= 0;
uh= 0;
um= 0;
uf=0;

a = 0;
i=0;
lamda =0;


v=0;


nf=0;
nfn=0;


y =0 ;
k = 0;
h = 0;
iv=0;
c = 0;
m = 0;
w =0 ;
r = 0;
rd = 0;
gg=0;
pf =0;
//q=0;
//eue=0;

yn =0 ;
kn = 0;
hn = 0;
cn = 0;
//mn = 0;
wn =0 ;
rn = 0;
rdn = 0;
ggn=0;
pfn =0;
pn = 0; 
pain = 0;
an = 0;
in=0;
ivn=0; 
lamdan=0;
//qn=0;


hg=0;
yg=0;

end;

@#endif


@#if tech==0

c-hb*c(-1)= (1+gamma*hb)*(c(+1)-hb*c) - gamma*hb*(c(+2)-hb*c(+1)) - ((1-hb)/sigc)*(1-gamma*hb)*((iss_noa/(1+iss_noa))*(1/iss_noa)*i + p - p(+1) + uc(+1) + ud(+1)) + ((1-hb)/sigc)*(uc + ud + gamma*hb*(uc(+2) +ud(+2))) ; 
sigh*h = w - uh + (1/(1-hb*gamma))*( (-sigc/(1-hb))*(c - hb*c(-1)) + uc ) + (hb*gamma/(1-hb*gamma))*( (sigc/(1-hb))*(c(+1)-hb*c) - uc(+1) - ud(+1) +ud );
sigm*(m-p) = (-1/(1+iss_noa))*(1/iss_noa)*i + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(c-hb*c(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(c(+1)-hb*c) + uc(+1) +ud(+1) - ud) );
lamda = (1-delta)*gamma*lamda(+1) - ((1-(1-delta)*gamma)/(1-hb*gamma))*( (sigc/(1-hb))*(c(+1)-hb*c) - (1/rss_noa)*r(+1) - uc(+1) - ud(+1) - hb*gamma*( (sigc/(1-hb))*(c(+2)-hb*c(+1)) - (1/rss_noa)*r(+1) -uc(+2) -ud(+2) ) );
iv =  (gamma/(1+gamma))*iv(+1) + (1/(1+gamma))*iv(-1) + (1/(ss*(1+gamma)))*lamda  + (1/(ss*(1+gamma)*(1-hb*gamma)))*( (sigc/(1-hb))*(c-hb*c(-1)) - uc - ud +gamma*hb*( (-sigc/(1-hb))*(c(+1)-hb*c) + uc(+1) +ud(+1) ) ) ;


yss_noa*y = css_noa*c + ivss_noa*iv +ggss_noa*gg;
(1-delta)*k(-1) +delta*iv - k = 0;

y = u + theta*k(-1) + (1-theta)*h;
h - k(-1) - (1/rss_noa)*r + w = 0;


@#if pindex==0
pai = gamma*(pai(+1)) + (((1-gamma*rho)*(1-rho))/(rho))*(-u + theta*(1/rss_noa)*r +(1-theta)*w ) ;
@#endif
@#if pindex==1
pai = ( gamma/(1+gamma*rho*index) )*(pai(+1)) + ( index/(1+gamma*rho*index) )*pai(-1)+ (  ((1-gamma*rho)*(1-rho))/( (rho)*(1+gamma*rho*index) )  )*(-u + theta*(1/rss_noa)*r +(1-theta)*w ); 
@#endif

pai = p - p(-1);


gg=(1/tau)*uf +y; // * govt expenditure

y_g=y-y(-1);
c_g=c-c(-1);
iv_g=iv-iv(-1);
h_g=h-h(-1);
w_g=w-w(-1);
tfp_g=tfp-tfp(-1);

tfpe= 0.999*tfpe(-1)+tfp_g + etfp;

tq=lamda-(1/(1-gamma*hb))*( (sigc/(1-hb))*(-c+hb*c(-1))+uc+ud - gamma*hb*( (sigc/(1-hb))*(-c(+1)+hb*c) + uc(+1)+ud(+1) ) );
sp=(tq+k(-1));
stk=sp;

//tq=-(lamda-(1/(1-gamma*hb))*( (sigc/(1-hb))*(-c+hb*c(-1))+uc+ud - gamma*hb*( (sigc/(1-hb))*(-c(+1)+hb*c) + uc(+1)+ud(+1) ) ));

//tq = sigc*c+lamda; //*Tobin's Q
//tq = -(css_noa^(-sigc))*sigc*c-(1/lamdass_noa)*lamda; //*Tobin's Q

rr=i-pai(+1); //*real interest rate
///////////////////////////////////////////////////////////////
// Monetary Policy (Inerest rate rule) & interest rate shock //
//////////////////////////////////////////////////////////////
iss_noa*(1/iss_noa)*i = ((1+iss_noa)^(xxx))*(xxx*(iss_noa/(1+iss_noa))*(1/iss_noa)*i(-1) + (1-xxx)*((xpai)*(pai)+(xy)*(yg)) + v);

@#if ishock==1
 v = rhoi*v(-1) + ei;
@#endif

@#if ishock==0
 v = rhoi*v(-1);
@#endif

@#if fshock==1
uf = rhouf*uf(-1)+ef;
@#endif

@#if fshock==0
uf = rhouf*uf(-1)+0.000001*ef;
@#endif

@#if ushock==1
 @#if news==1
   
  
   @#if eshock==16
   u = rhou*u(-1) +eu + eue1(-1) + eue2(-2) + eue3(-3) + eue4(-4)+ eue5(-5)+ eue6(-6)+ eue7(-7)+ eue8(-8)+ eue9(-9) +eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue(-16); 
   
   
   @#if var_data==22  
    conf12= eue(-12); conf11= eue11(-11);conf10= eue10(-10);conf9= eue9(-9);conf8= eue8(-8);
    //conf=rhou*(rhou^1*(u)+eue1)+eue2;
    conf=((rhou^1*(u)+eue1)
        + rhou*(rhou^1*(u)+eue1)+eue2  
        + rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3
        +  rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4
        +  rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9 
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10 
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11)+eue12
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11)+eue12)+eue13
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11)+eue12)+eue13)+eue14
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11)+eue12)+eue13)+eue14)+eue15
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11)+eue12)+eue13)+eue14)+eue15)+eue 
        )/1;
   @#endif

   
   @#endif
   
 @#endif
 @#if news==0
   

 @#if eshock==16
    u = rhou*u(-1); 
     @#if var_data==20
    conf12= 0; conf11=0;conf10= 0;conf9=0;conf8= 0;
    @#if conf_var==1
    conf=0;
   @#endif
   @#if conf_var==2
      conf=0;
   @#endif
    @#endif

   @#if var_data==21
    conf12= 0; conf11=0;conf10= 0;conf9=0;conf8= 0;
   @#if conf_var==1
     conf=0;
     @#endif
     @#if conf_var==2
       conf=0;
     @#endif
    @#endif


  @#if var_data==25
    conf12= 0; conf11=0;conf10= 0;conf9=0;conf8= 0;
    @#if conf_var==1
    conf=0;
   @#endif
   @#if conf_var==2
      conf=0;
   @#endif
    @#endif
    
    @#if var_data==22 
    conf12= 0; conf11=0;conf10= 0;conf9=0;conf8= 0;
    @#if conf_var==1
    conf=0;
   @#endif
   @#if conf_var==2
    conf=0;
   @#endif

  @#endif 
   @#endif  
 



 @#endif
@#endif 

@#if ushock==0
 @#if news==1
   
   @#if eshock==16
   u = rhou*u(-1) + eue1(-1) + eue2(-2) + eue3(-3) + eue4(-4)+ eue5(-5)+ eue6(-6)+ eue7(-7)+ eue8(-8)+ eue9(-9) +eue10(-10)+eue11(-11)+eue12(-12)+eue13(-13)+eue14(-14)+eue15(-15)+eue(-16); 
   
   @#if var_data==22  
    conf12= eue(-12); conf11= eue11(-11);conf10= eue10(-10);conf9= eue9(-9);conf8= eue8(-8);
    conf=((rhou^1*(u)+eue1)
        + rhou*(rhou^1*(u)+eue1)+eue2  
        + rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3
        +  rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4
        +  rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9 
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10 
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11)+eue12
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11)+eue12)+eue13
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11)+eue12)+eue13)+eue14
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11)+eue12)+eue13)+eue14)+eue15
        +  rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou*(rhou^1*(u)+eue1)+eue2)+eue3)+eue4)+eue5)+eue6)+eue7)+eue8)+eue9)+eue10)+eue11)+eue12)+eue13)+eue14)+eue15)+eue 
        )/1;
   @#endif

  @#endif

 @#endif


 @#if news==0



 @#if eshock==16
    u = rhou*u(-1); 
       
    @#if var_data==22 
    conf12= 0; conf11=0;conf10= 0;conf9=0;conf8= 0;
    @#if conf_var==1
    conf=0;
   @#endif
   @#if conf_var==2
    conf=0;
   @#endif

  @#endif 
   @#endif  
  


 @#endif 
@#endif

 
////////////////////////
// Preference shocks //
//////////////////////

ud = rhoud*ud(-1); //* Preference shock


@#if cshock==1
uc = rhouc*uc(-1) +euc; //* Preference shock
@#endif

@#if cshock==0
uc = rhouc*uc(-1); //* Preference shock
@#endif


@#if hshock==1
uh = rhouh*uh(-1) +euh; //* Preference shock
@#endif

@#if hshock==0
uh = rhouh*uh(-1); //* Preference shock
@#endif


//um = rhoum*um(-1) +eum; //* Preference shock
um = rhoum*um(-1) ; //* Preference shock


//yg=y-yn; //* gdp gap
yg=y-ynobs; //* gdp gap with measurement errors
hg = h-hn;


mcgap=(-u + theta*(1/rss_noa)*r +(1-theta)*w );
mcgape=mcgap+mee;
mc=mcgap;
mce=mcgape;
//gap_nw=(-u + theta*(1/rss_noa)*r +(1-theta)*w );
fpai = pai(+1);
ppai = pai(-1);
fpaie = pai(+1)+meee;
paiee=pai+meee+me;//*pai with 2sls errors & measument errors
paie=pai+me;//*pai with measument errors
ppaie=paie(-1);

tfp=u;
tfpu=u;
shu=eu;
shi=ei;
shc=euc;
shh=euh;
shue1=eue1;

/////////////////////////
//* natural variables *//
////////////////////////


cn-hb*cn(-1)= (1+gamma*hb)*(cn(+1)-hb*cn) - gamma*hb*(cn(+2)-hb*cn(+1)) - ((1-hb)/sigc)*(1-gamma*hb)*((iss_noa/(1+iss_noa))*(1/iss_noa)*in + pn - pn(+1) + uc(+1) + ud(+1)) + ((1-hb)/sigc)*(uc + ud + gamma*hb*(uc(+2) +ud(+2))) ; 
sigh*hn = wn - uh + (1/(1-hb*gamma))*( (-sigc/(1-hb))*(cn - hb*cn(-1)) + uc ) + (hb*gamma/(1-hb*gamma))*( (sigc/(1-hb))*(cn(+1)-hb*cn) - uc(+1) - ud(+1) +ud );
//sigm*(mn-pn) = (-1/(1+iss_noa))*(1/iss_noa)*in + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(cn-hb*cn(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) - ud) );
sigm*(0-pn) = (-1/(1+iss_noa))*(1/iss_noa)*in + um + (1/(1-hb*gamma))*( (sigc/(1-hb))*(cn-hb*cn(-1)) -uc + (hb*gamma)*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) - ud) );

lamdan = (1-delta)*gamma*lamdan(+1) - ((1-(1-delta)*gamma)/(1-hb*gamma))*( (sigc/(1-hb))*(cn(+1)-hb*cn) - (1/rss_noa)*rn(+1) - uc(+1) - ud(+1) - hb*gamma*( (sigc/(1-hb))*(cn(+2)-hb*cn(+1)) - (1/rss_noa)*rn(+1) -uc(+2) -ud(+2) ) );
ivn =  (gamma/(1+gamma))*ivn(+1) + (1/(1+gamma))*ivn(-1) + (1/(ss*(1+gamma)))*lamdan  + (1/(ss*(1+gamma)*(1-hb*gamma)))*( (sigc/(1-hb))*(cn-hb*cn(-1)) - uc - ud +gamma*hb*( (-sigc/(1-hb))*(cn(+1)-hb*cn) + uc(+1) +ud(+1) ) ) ;

 yss_noa*yn = css_noa*cn + ivss_noa*ivn +ggss_noa*ggn;
(1-delta)*kn(-1) +delta*ivn - kn = 0;
//ynobs=yn+me;
ynobs=yn;

hn - kn(-1) - (1/rss_noa)*rn + wn = 0;
yn = u + theta*kn(-1) + (1-theta)*hn;
0 = (((1-gamma*0)*(1-0)))*(-u + theta*(1/rss_noa)*rn +(1-theta)*wn ) ;

pain = pn - pn(-1);
//mn = mn(-1);
//mn = 0;

//tqn = sigc*cn+lamdan; //*Tobin's Q
tfpn=u;
tqn=lamdan-(1/(1-gamma*hb))*( (sigc/(1-hb))*(-cn+hb*cn(-1))+uc+ud - gamma*hb*( (sigc/(1-hb))*(-cn(+1)+hb*cn) + uc(+1)+ud(+1) ) );
//tqn=-(lamdan-(1/(1-gamma*hb))*( (sigc/(1-hb))*(-cn+hb*cn(-1))+uc+ud - gamma*hb*( (sigc/(1-hb))*(-cn(+1)+hb*cn) + uc(+1)+ud(+1) ) ));


mcn=(-u + theta*(1/rss_noa)*rn +(1-theta)*wn );

ggn= (1/tau)*uf +yn;

end;


initval;

//sh=0;
p = 0; 
pai = 0;
u = 0;
ud= 0;
uc= 0;
uh= 0;
um= 0;


i=0;
lamda =0;
iv=0;


v=0;



y =0 ;
k = 0;
h = 0;
iv=0;
c = 0;
m = 0;
w =0 ;
r = 0;

//q=0;
//eue=0;

yn =0 ;
kn = 0;
hn = 0;
cn = 0;
//mn = 0;
wn =0 ;
rn = 0;
//rdn = 0;

pn = 0; 
pain = 0;

in=0;
ivn=0; 
lamdan=0;
//qn=0;


hg=0;
yg=0;

end;

steady_state_model;
p = 0; 
pai = 0;
u = 0;
ud= 0;
uc= 0;
uh= 0;
um= 0;


i=0;
lamda =0;
//tiv=0;

v=0;




y =0 ;
k = 0;
h = 0;
iv=0;
c = 0;
m = 0;
w =0 ;
r = 0;

//q=0;
//eue=0;

yn =0 ;
kn = 0;
hn = 0;
cn = 0;
//mn = 0;
wn =0 ;
rn = 0;

//pfn =0;
pn = 0; 
pain = 0;

in=0;
ivn=0; 
lamdan=0;
//qn=0;


hg=0;
yg=0;
end;
@#endif


 shocks;

 var me;  // measurement errors  for 'mc'
 stderr me_std; 

 var mee;  // measurement errors  for 'mc'
 stderr mee_std;

 var meee;  // measurement errors for 'fpai'
 stderr meee_std;




//***********************************//
//**** Start of fig==0 Section **////
//***********************************//
@#if fig==0
@#if data_techzero==0  
 var eu;
 stderr aa0(j,1); 
 
 var ei;
 stderr aa0(j,2);

 var euc;
 stderr aa0(j,3);

  var euh;
 stderr aa0(j,4);

  var ef;
 stderr aa0(j,5);

 var etfp;
 stderr tfpwt*aa0(j,10);

@#endif

@#if data_techzero==1  
 var eu;
stderr aa0(j,1); 
//stderr 0; 
 
 var ei;
@#if var_shock==0
stderr aa0(j,2);
@#endif
@#if var_shock==1
stderr 0; 
@#endif
 
var euc;
@#if var_shock==0
stderr aa0(j,3);
@#endif
@#if var_shock==1
stderr 0; 
@#endif 

var euh;
@#if var_shock==0
stderr aa0(j,4);
@#endif
@#if var_shock==1
stderr 0; 
@#endif 


var ef;
@#if var_shock==0
stderr aa0(j,5);
@#endif
@#if var_shock==1
stderr 0; 
@#endif 

var etfp;
 stderr tfpwt*aa0(j,10);

@#endif
///////////////////
// NEWS SHOCKS//
//////////////////
@#if data_techzero==0 

  var eue1;
 stderr 0;

  var eue2;
 stderr 0;

  var eue3;
 stderr 0;

  var eue4;
 stderr aa0(j,6);

  var eue5;
 stderr 0;

  var eue6;
 stderr 0;

  var eue7;
 stderr  0;

 var eue8;
 stderr aa0(j,7);

  var eue9;
 stderr  0;
 
  var eue10;
 stderr  0;
   

@#if eshock==16 
 var eue11;
 stderr 0;

  var eue12;
 stderr aa0(j,8);

 var eue13;
 stderr 0;
 
 var eue14;
 stderr 0;

 var eue15;
 stderr 0;
 
 var eue;
 stderr aa0(j,9);


@#if crosscrr==1
corr eue4, eue8 =aa0(j,11); corr eue4, eue12 =aa0(j,12); corr eue4, eue =aa0(j,13); corr eue4, eu =0;
corr eue8, eue12= aa0(j,14); corr eue8, eue =aa0(j,15); corr eue8, eu =0;
corr eue12, eue =aa0(j,16); corr eue12, eu =0;
corr eue, eu =0;
@#endif



@#endif

   
@#endif

@#if data_techzero==1 

  var eue1;
 stderr 0;

  var eue2;
 stderr 0;

  var eue3;
 stderr 0;

  var eue4;
 stderr aa0(j,6);

  var eue5;
 stderr 0;

  var eue6;
 stderr 0;

  var eue7;
 stderr 0;

 var eue8;
 stderr aa0(j,7);

  var eue9;
 stderr 0;

  var eue10;
  stderr 0;



@#if eshock==16
 var eue11;
 stderr 0;

  var eue12;
 stderr aa0(j,8);

 var eue13;
 stderr 0;
 
 var eue14;
 stderr 0;

 var eue15;
 stderr 0;
 
 var eue;
 stderr aa0(j,9);


@#if crosscrr==1
corr eue4, eue8 =aa0(j,11); corr eue4, eue12 =aa0(j,12); corr eue4, eue =aa0(j,13); corr eue4, eu =0;
corr eue8, eue12 = aa0(j,14); corr eue8, eue =aa0(j,15); corr eue8, eu =0;
corr eue12, eue =aa0(j,16); corr eue12, eu =0;

//corr eue4, eue8 =0.9; corr eue4, eue12 =0.411; corr eue4, eue =0.18; corr eue4, eu =0;
//corr eue8, eue12 = 0.445; corr eue8, eue =0.188; corr eue8, eu =0;
//corr eue12, eue =0.95; corr eue12, eu =0;

corr eue, eu =0;
@#endif

@#endif

 

@#endif

  
@#endif

//@#endif
//***********************************//
//** End of fig==0 section **//
//***********************************//


 
end;



for ii=iia:iib
set_dynare_seed(ii); 


@#if large_sample==0

@#if tech==1
 @#if fig==0
  @#if hp==0
    
    @#if var_data==22
    stoch_simul(irf=40, order=1,periods=6000,nograph, noprint ,qz_criterium=0.999999999, drop=100)  y a u c iv tiv k h rd pai fpai p tq mc nf  pf_exp w r rr i tfp tfpe tfpu  yn an cn ivn tivn kn hn rdn pain tqn mcn nfn sp conf pf stk y_g  c_g  tiv_g iv_g  h_g  w_g tfp_g;    
     @#endif
    @#endif
 @#endif
 
 @#if fig==1
 @#if cndvar==0
  @#if full_fig==0 
  stoch_simul(irf=40, order=1,periods=580,nograph, drop=10) y a u c iv tiv k h rd pai fpai p  tq mc nf  pf pf_exp w r rr i tfp tfpu  yn an cn ivn tivn kn hn rdn pain tqn mcn nfn sp  conf pf stk y_g  c_g  tiv_g iv_g  h_g  w_g tfp_g;  
  @#endif
  @#if full_fig==1 
  stoch_simul(irf=40, order=1,periods=580, drop=10) y a u c iv tiv k h rd pai fpai p  tq mc nf  pf pf_exp w r rr i tfp tfpu  yn an cn ivn tivn kn hn rdn pain tqn mcn nfn sp  conf pf stk y_g  c_g  tiv_g iv_g  h_g  w_g tfp_g;  
  @#endif 
 @#endif 
 @#if cndvar==1
 stoch_simul(irf=40, order=1,periods=0,nograph, drop=10,conditional_variance_decomposition = [ 8 12 16 20 40]) y a u c iv tiv k h rd pai fpai p  tq mc nf   pf pf_exp w r rr i tfp tfpu  yn an cn ivn tivn kn hn rdn pain tqn mcn nfn sp  conf pf stk y_g  c_g  tiv_g iv_g  h_g  w_g tfp_g;  
 @#endif
@#endif
@#endif

@#if tech==0
 @#if fig==0
   @#if hp==0
    @#if var_data==22
    stoch_simul(irf=40, order=2,periods=6000,nograph,noprint ,qz_criterium=0.999999, drop=100)  y  u mc c iv  k  h  pai fpai p  w r rr i tfp tfpe  tfpu  tq   mcn sp  conf   stk y_g  c_g  h_g iv_g  w_g tfp_g; 
    @#endif
        
  @#endif
  @#if hp==1
   stoch_simul(irf=40, order=1,periods=580,nograph,hp_filter=1600, drop=10)   y yn u mc c iv  k  h  pai fpai  p w r  rr i tfp tfpu  tq tqn yn cn ivn kn hn pain tqn mcn sp  conf   stk; 
  @#endif
 @#endif   
 @#if fig==1
 @#if cndvar==0
  @#if full_fig==0 
  stoch_simul(irf=40, order=1,periods=580,nograph, drop=10)  y yn u mc c iv  k  h  pai fpai p  w r rr i tfp tfpu  tq tqn yn cn ivn kn hn pain tqn mcn sp  conf  stk; 
  @#endif
  @#if full_fig==1  
  stoch_simul(irf=40, order=1,periods=580, drop=10)  y yn u mc c iv  k  h  pai fpai p  w r rr i tfp tfpu  tq tqn yn cn ivn kn hn pain tqn mcn sp  conf stk; 
  @#endif
 @#endif
 @#if cndvar==1
 stoch_simul(irf=40, order=1,periods=0,nograph, drop=10,conditional_variance_decomposition = [1 4 8 12 16 20 40])  y yn u mc c iv  k  h  pai fpai p  w r rr i tfp tfpu  tq tqn yn cn ivn kn hn pain tqn mcn  sp  conf  stk; 
 @#endif
 @#endif   
@#endif

@#endif



//*******************************************************//
//*** Creating simulated dataset for monte carolo VAR ***// 
//*******************************************************//

@#if sim_var==1
 
/* disricarding first 400 obs*/
@#if tech==1

@#if bps==0

@#if trend_sigma==0
 @#if tfp_errors==1
 var_monte_data1{:,1}=tfpe(inds:inde);
 @#endif
 @#if tfp_errors==0
 var_monte_data1{:,1}=tfp(inds:inde);
 @#endif
var_monte_data1{:,2}=y(inds:inde);
var_monte_data1{:,3}=c(inds:inde);
var_monte_data1{:,4}=w(inds:inde);
//var_monte_data1{:,4}=iv(2101:2291);
var_monte_data1{:,5}=h(inds:inde);
@#endif


var_monte_data1{:,6}=i(inds:inde);
var_monte_data1{:,7}=pai(inds:inde);


@#if var_data==21
 @#if eshock==16
 var_monte_data1{:,11}=(conf(2000:2190));
 @#endif
@#endif


@#endif


@#endif

@#if tech==0

@#if bps==0

@#if trend_sigma==0
@#if tfp_errors==1
var_monte_data1{:,1}=tfpe(inds:inde);
@#endif
@#if tfp_errors==0
var_monte_data1{:,1}=tfp(inds:inde);
@#endif
var_monte_data1{:,2}=y(inds:inde);
var_monte_data1{:,3}=c(inds:inde);
//var_monte_data1{:,4}=tiv(inds:inde);
var_monte_data1{:,4}=w(inds:inde);
var_monte_data1{:,5}=h(inds:inde);
@#endif


var_monte_data1{:,6}=i(inds:inde);
var_monte_data1{:,7}=pai(inds:inde);

 @#if var_data==21
 var_monte_data1{:,11}=conf(inds:inde);
 @#endif


@#endif

@#endif

var_monte_data=cell2mat(var_monte_data1);
var_monte_data_i{:,:,ii}=var_monte_data;
var_monte_data_ii=cell2mat(var_monte_data_i);
save 'var_monte_data' var_monte_data;
%save 'var_monte_data_ii'

KurmannOtrok_main_sw_var3_cp2;

@#if var_est_method==2
 
 @#if var_data==22
var_monte_tfp_med(:,ii,j)=datamed(:,7);
var_monte_y_med(:,ii,j)=datamed(:,8);
var_monte_c_med(:,ii,j)=datamed(:,9);
var_monte_h_med(:,ii,j)=datamed(:,10);
var_monte_i_med(:,ii,j)=datamed(:,11);
var_monte_pai_med(:,ii,j)=datamed(:,12);
@#endif 

 
 
@#endif


@#if var_est_method==1
 

@#if var_data==22
var_monte_tfp_med(:,ii,j)=datamed(:,7);
var_monte_y_med(:,ii,j)=datamed(:,8);
var_monte_c_med(:,ii,j)=datamed(:,9);
var_monte_h_med(:,ii,j)=datamed(:,10);
var_monte_i_med(:,ii,j)=datamed(:,11);
var_monte_pai_med(:,ii,j)=datamed(:,12);
@#endif


@#endif

@#if var_est_method==0
 

@#if var_data==22
var_monte_tfp_med(:,ii,j)=data_k(:,7);
var_monte_y_med(:,ii,j)=data_k(:,8);
var_monte_c_med(:,ii,j)=data_k(:,9);
var_monte_h_med(:,ii,j)=data_k(:,10);
var_monte_i_med(:,ii,j)=data_k(:,11);
var_monte_pai_med(:,ii,j)=data_k(:,12);
@#endif

  
@#endif

@#endif
//********************************************************//


end
//** End of 'for ii' : 'ii' is for a monte carlo sim, i.e., shocks are drawn "ii" times **//



end 
//**************************************************************************************************************************************//
//** End of 'for j' : 'j' is for Bayesian coefficient estimates, i.e., different coeffifient estimates (j different coeffs) are used **//
//**************************************************************************************************************************************//



//*******************************************************************************************************//
//*** Start of Monte carlo (bayesian?) multiple news shocks IRF (can use Bayesian by choosing e.g. ja=1 &jb=200 *****//
//*******************************************************************************************************//




@#if sim_var==1

@#if var_est_method==1


@#if var_data==22

size_mat=size(var_monte_tfp_med);
rw=size_mat(1);
cl=size_mat(2)*size_mat(3);

var_monte_tfp_med_rsp=reshape(var_monte_tfp_med,[rw cl]);
var_monte_y_med_rsp=reshape(var_monte_y_med,[rw cl]);
var_monte_c_med_rsp=reshape(var_monte_c_med,[rw cl]);
var_monte_h_med_rsp=reshape(var_monte_h_med,[rw cl]);
var_monte_i_med_rsp=reshape(var_monte_i_med,[rw cl]);
var_monte_pai_med_rsp=reshape(var_monte_pai_med,[rw cl]);

@#if median_mean==1
var_monte_tfp_av_av=median(var_monte_tfp_med_rsp,2);
var_monte_y_av_av=median(var_monte_y_med_rsp,2);
var_monte_c_av_av=median(var_monte_c_med_rsp,2);
var_monte_h_av_av=median(var_monte_h_med_rsp,2);
var_monte_i_av_av=median(var_monte_i_med_rsp,2);
var_monte_pai_av_av=median(var_monte_pai_med_rsp,2);
@#endif

@#if median_mean==2
var_monte_tfp_av_av=mean(var_monte_tfp_med_rsp,2);
var_monte_y_av_av=mean(var_monte_y_med_rsp,2);
var_monte_c_av_av=mean(var_monte_c_med_rsp,2);
var_monte_h_av_av=mean(var_monte_h_med_rsp,2);
var_monte_i_av_av=mean(var_monte_i_med_rsp,2);
var_monte_pai_av_av=mean(var_monte_pai_med_rsp,2);
@#endif


var_monte_tfp_pt5=prctile(var_monte_tfp_med_rsp,5,2);
var_monte_tfp_pt95=prctile(var_monte_tfp_med_rsp,95,2);
var_monte_y_pt5=prctile(var_monte_y_med_rsp,5,2);
var_monte_y_pt95=prctile(var_monte_y_med_rsp,95,2);
var_monte_c_pt5=prctile(var_monte_c_med_rsp,5,2);
var_monte_c_pt95=prctile(var_monte_c_med_rsp,95,2);
var_monte_h_pt5=prctile(var_monte_h_med_rsp,5,2);
var_monte_h_pt95=prctile(var_monte_h_med_rsp,95,2);
var_monte_i_pt5=prctile(var_monte_i_med_rsp,5,2);
var_monte_i_pt95=prctile(var_monte_i_med_rsp,95,2);
var_monte_pai_pt5=prctile(var_monte_pai_med_rsp,5,2);
var_monte_pai_pt95=prctile(var_monte_pai_med_rsp,95,2);

var_monte_tfp_ptstdl=prctile(var_monte_tfp_med_rsp,16.5,2);
var_monte_tfp_ptstdh=prctile(var_monte_tfp_med_rsp,83.5,2);
var_monte_y_ptstdl=prctile(var_monte_y_med_rsp,16.5,2);
var_monte_y_ptstdh=prctile(var_monte_y_med_rsp,83.5,2);
var_monte_c_ptstdl=prctile(var_monte_c_med_rsp,16.5,2);
var_monte_c_ptstdh=prctile(var_monte_c_med_rsp,83.5,2);
var_monte_h_ptstdl=prctile(var_monte_h_med_rsp,16.5,2);
var_monte_h_ptstdh=prctile(var_monte_h_med_rsp,83.5,2);
var_monte_i_ptstdl=prctile(var_monte_i_med_rsp,16.5,2);
var_monte_i_ptstdh=prctile(var_monte_i_med_rsp,83.5,2);
var_monte_pai_ptstdl=prctile(var_monte_pai_med_rsp,16.5,2);
var_monte_pai_ptstdh=prctile(var_monte_pai_med_rsp,83.5,2);






@#endif




@#endif

@#if var_est_method==2

@#if var_data==21
var_monte_tfp_av_av=median(var_monte_tfp_med,2);
var_monte_y_av_av=median(var_monte_y_med,2);
var_monte_c_av_av=median(var_monte_c_med,2);
var_monte_h_av_av=median(var_monte_h_med,2);
var_monte_i_av_av=median(var_monte_i_med,2);
var_monte_pai_av_av=median(var_monte_pai_med,2);
var_monte_conf_av_av=median(var_monte_conf_med,2);
 @#endif

@#if var_data==22
//var_monte_tfp_av_av=median(var_monte_tfp_med,2);
//var_monte_y_av_av=median(var_monte_y_med,2);
//var_monte_c_av_av=median(var_monte_c_med,2);
//var_monte_h_av_av=median(var_monte_h_med,2);
//var_monte_i_av_av=median(var_monte_i_med,2);
//var_monte_pai_av_av=median(var_monte_pai_med,2);




size_mat=size(var_monte_tfp_med);
rw=size_mat(1);
cl=size_mat(2)*size_mat(3);


var_monte_tfp_med_rsp=reshape(var_monte_tfp_med,[rw cl]);
var_monte_y_med_rsp=reshape(var_monte_y_med,[rw cl]);
var_monte_c_med_rsp=reshape(var_monte_c_med,[rw cl]);
var_monte_h_med_rsp=reshape(var_monte_h_med,[rw cl]);
var_monte_i_med_rsp=reshape(var_monte_i_med,[rw cl]);
var_monte_pai_med_rsp=reshape(var_monte_pai_med,[rw cl]);

@#if median_mean==1
var_monte_tfp_av_av=median(var_monte_tfp_med_rsp,2);
var_monte_y_av_av=median(var_monte_y_med_rsp,2);
var_monte_c_av_av=median(var_monte_c_med_rsp,2);
var_monte_h_av_av=median(var_monte_h_med_rsp,2);
var_monte_i_av_av=median(var_monte_i_med_rsp,2);
var_monte_pai_av_av=median(var_monte_pai_med_rsp,2);
@#endif

@#if median_mean==2
var_monte_tfp_av_av=mean(var_monte_tfp_med_rsp,2);
var_monte_y_av_av=mean(var_monte_y_med_rsp,2);
var_monte_c_av_av=mean(var_monte_c_med_rsp,2);
var_monte_h_av_av=mean(var_monte_h_med_rsp,2);
var_monte_i_av_av=mean(var_monte_i_med_rsp,2);
var_monte_pai_av_av=mean(var_monte_pai_med_rsp,2);
@#endif

var_monte_tfp_pt5=prctile(var_monte_tfp_med_rsp,5,2);
var_monte_tfp_pt95=prctile(var_monte_tfp_med_rsp,95,2);
var_monte_y_pt5=prctile(var_monte_y_med_rsp,5,2);
var_monte_y_pt95=prctile(var_monte_y_med_rsp,95,2);
var_monte_c_pt5=prctile(var_monte_c_med_rsp,5,2);
var_monte_c_pt95=prctile(var_monte_c_med_rsp,95,2);
var_monte_h_pt5=prctile(var_monte_h_med_rsp,5,2);
var_monte_h_pt95=prctile(var_monte_h_med_rsp,95,2);
var_monte_i_pt5=prctile(var_monte_i_med_rsp,5,2);
var_monte_i_pt95=prctile(var_monte_i_med_rsp,95,2);
var_monte_pai_pt5=prctile(var_monte_pai_med_rsp,5,2);
var_monte_pai_pt95=prctile(var_monte_pai_med_rsp,95,2);

var_monte_tfp_ptstdl=prctile(var_monte_tfp_med_rsp,16.5,2);
var_monte_tfp_ptstdh=prctile(var_monte_tfp_med_rsp,83.5,2);
var_monte_y_ptstdl=prctile(var_monte_y_med_rsp,16.5,2);
var_monte_y_ptstdh=prctile(var_monte_y_med_rsp,83.5,2);
var_monte_c_ptstdl=prctile(var_monte_c_med_rsp,16.5,2);
var_monte_c_ptstdh=prctile(var_monte_c_med_rsp,83.5,2);
var_monte_h_ptstdl=prctile(var_monte_h_med_rsp,16.5,2);
var_monte_h_ptstdh=prctile(var_monte_h_med_rsp,83.5,2);
var_monte_i_ptstdl=prctile(var_monte_i_med_rsp,16.5,2);
var_monte_i_ptstdh=prctile(var_monte_i_med_rsp,83.5,2);
var_monte_pai_ptstdl=prctile(var_monte_pai_med_rsp,16.5,2);
var_monte_pai_ptstdh=prctile(var_monte_pai_med_rsp,83.5,2);

@#endif




@#endif

@#if var_est_method==0


@#if var_data==22


size_mat=size(var_monte_tfp_med);
rw=size_mat(1);
cl=size_mat(2)*size_mat(3);

var_monte_tfp_med_rsp=reshape(var_monte_tfp_med,[rw cl]);
var_monte_y_med_rsp=reshape(var_monte_y_med,[rw cl]);
var_monte_c_med_rsp=reshape(var_monte_c_med,[rw cl]);
var_monte_h_med_rsp=reshape(var_monte_h_med,[rw cl]);
var_monte_i_med_rsp=reshape(var_monte_i_med,[rw cl]);
var_monte_pai_med_rsp=reshape(var_monte_pai_med,[rw cl]);


@#if median_mean==1
var_monte_tfp_av_av=median(var_monte_tfp_med_rsp,2);
var_monte_y_av_av=median(var_monte_y_med_rsp,2);
var_monte_c_av_av=median(var_monte_c_med_rsp,2);
var_monte_h_av_av=median(var_monte_h_med_rsp,2);
var_monte_i_av_av=median(var_monte_i_med_rsp,2);
var_monte_pai_av_av=median(var_monte_pai_med_rsp,2);
@#endif

@#if median_mean==2
var_monte_tfp_av_av=mean(var_monte_tfp_med_rsp,2);
var_monte_y_av_av=mean(var_monte_y_med_rsp,2);
var_monte_c_av_av=mean(var_monte_c_med_rsp,2);
var_monte_h_av_av=mean(var_monte_h_med_rsp,2);
var_monte_i_av_av=mean(var_monte_i_med_rsp,2);
var_monte_pai_av_av=mean(var_monte_pai_med_rsp,2);
@#endif



var_monte_tfp_pt5=prctile(var_monte_tfp_med_rsp,5,2);
var_monte_tfp_pt95=prctile(var_monte_tfp_med_rsp,95,2);
var_monte_y_pt5=prctile(var_monte_y_med_rsp,5,2);
var_monte_y_pt95=prctile(var_monte_y_med_rsp,95,2);
var_monte_c_pt5=prctile(var_monte_c_med_rsp,5,2);
var_monte_c_pt95=prctile(var_monte_c_med_rsp,95,2);
var_monte_h_pt5=prctile(var_monte_h_med_rsp,5,2);
var_monte_h_pt95=prctile(var_monte_h_med_rsp,95,2);
var_monte_i_pt5=prctile(var_monte_i_med_rsp,5,2);
var_monte_i_pt95=prctile(var_monte_i_med_rsp,95,2);
var_monte_pai_pt5=prctile(var_monte_pai_med_rsp,5,2);
var_monte_pai_pt95=prctile(var_monte_pai_med_rsp,95,2);

var_monte_tfp_ptstdl=prctile(var_monte_tfp_med_rsp,16.5,2);
var_monte_tfp_ptstdh=prctile(var_monte_tfp_med_rsp,83.5,2);
var_monte_y_ptstdl=prctile(var_monte_y_med_rsp,16.5,2);
var_monte_y_ptstdh=prctile(var_monte_y_med_rsp,83.5,2);
var_monte_c_ptstdl=prctile(var_monte_c_med_rsp,16.5,2);
var_monte_c_ptstdh=prctile(var_monte_c_med_rsp,83.5,2);
var_monte_h_ptstdl=prctile(var_monte_h_med_rsp,16.5,2);
var_monte_h_ptstdh=prctile(var_monte_h_med_rsp,83.5,2);
var_monte_i_ptstdl=prctile(var_monte_i_med_rsp,16.5,2);
var_monte_i_ptstdh=prctile(var_monte_i_med_rsp,83.5,2);
var_monte_pai_ptstdl=prctile(var_monte_pai_med_rsp,16.5,2);
var_monte_pai_ptstdh=prctile(var_monte_pai_med_rsp,83.5,2);

@#endif



@#endif
@#endif

@#if fignum==400
figure(400)
@#endif
@#if fignum==500
figure(500)
@#endif


@#if sim_var==1

subplot(3,3,1)
if lsp ==1
plot(var_monte_tfp_av_av,'r--');
elseif lsp ==2
plot(var_monte_tfp_av_av,'r.-');
end
title('tfp');
hold on

 
subplot(3,3,2)
if lsp ==1
plot(var_monte_y_av_av,'r--');
elseif lsp ==2
plot(var_monte_y_av_av,'r.-');
end
title('y');
hold on


subplot(3,3,3)
if lsp ==1
plot(var_monte_c_av_av,'r--');
elseif lsp ==2
plot(var_monte_c_av_av,'r.-');
end
title('c');
hold on

 
@#if var_data==22
subplot(3,3,5)
if lsp ==1
plot(var_monte_h_av_av,'r--');
elseif lsp ==2
plot(var_monte_h_av_av,'r.-');
end
title('h');
hold on
@#endif

@#if var_data==22
subplot(3,3,6)
if lsp ==1
plot(var_monte_i_av_av,'r--');
elseif lsp ==2
plot(var_monte_i_av_av,'r.-');
end
title('i');
hold on
@#endif


subplot(3,3,7)
if lsp ==1
plot(var_monte_pai_av_av,'r--');
elseif lsp ==2
plot(var_monte_pai_av_av,'r.-');
end
title('pai');
hold on


@#if var_data==0
subplot(3,3,8)
if lsp ==1
plot(var_monte_sp_av_av,'r--');
elseif lsp ==2
plot(var_monte_sp_av_av,'r.-');
end
title('stock');
hold on
 @#endif



@#endif







